FAST COMPUTATION OF HIGH FREQUENCY DIRICHLET 
EIGENMODES VIA THE SPECTRAL FLOW OF THE INTERIOR 
NEUMANN-TO-DIRICHLET MAP 



ALEX BARNETT AND ANDREW HASSELL 



Abstract. We present a new algorithm for numerical computation of large 
eigenvalues and associated eigenfunctions of the Dirichlet Laplacian in a smooth, 
star-shaped domain in R d , d > 2. Conventional boundary-based methods re- 
quire a root-search in eigenfrcquency k, hence take 0(7V 3 ) effort per eigen- 
pair found, using dense linear algebra, where TV = 0(k d ~ 1 ) is the number of 
unknowns required to discretize the boundary. Our method is O(N) faster, 
achieved by linearizing with respect to k the spectrum of a weighted interior 
Neumann-to-Dirichlet (NtD) operator for the Hclmholtz equation. Approxi- 
mations kj to the square-roots kj of all O(N) eigenvalues lying in [k — e,k], 
where e = O(l), are found with 0(N S ) effort. We prove an error estimate 

l%-fcil<c(^ + £ 3 ), 

with C independent of k. We present a higher-order variant with eigenvalue 
error scaling empirically as 0(e 5 ) and eigenfunction error as 0(e 3 ), the former 
improving upon the 'scaling method' of Vergini— Saraceno. For planar domains 
(d = 2), with an assumption of absence of spectral concentration, we also prove 
rigorous error bounds that are close to those numerically observed. For d = 2 
we compute robustly the spectrum of the NtD operator via potential theory, 
Nystrom discretization, and the Cayley transform. At high frequencies (400 
wavelengths across), with eigenfrequency relative error 10 — 10 , we show that 
the method is 10 3 times faster than standard ones based upon a root-search. 
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1. Introduction 

Let f2 be a smooth, bounded domain in R d , strictly star-shaped with respect to 
the origin, that is x ■ n > for each x € dQ where n is the outward-pointing unit 
normal vector. We are interested in computing numerically the eigenvalues kj, and 
eigenfunctions or eigenmodes <j)j (normalized by H^Hl 2 ^) — 1)) of the Dirichlet 
Laplacian A = J2Li 5 V^ s on O. That is, 

(A + = in n , (1) 

4> =0 onffl. (2) 

We will refer to kj, the square-roots of eigenvalues, as (Dirichlet) cigcnfrequencies, 
and order them < fci < fc 2 < k 3 < . . . counting multiplicities. This classical 
problem has many applications in engineering and physics )20lH], principally in the 
modeling of acoustic, electromagnetic and optical cavities, vibrating membranes, 
trapped quantum particles and nano-scale devices [44], and in data analysis |51) . 
Note that some applications involve homogeneous boundary conditions other than 
([2]), or the Maxwell or elasticity equations, yet the above serves as a paradigm for 
this larger class of problems. In d = 2 it is known as the 'drum' problem, and is 
reviewed in |39[ 156] . A numerical approach is needed for all but the small subset 
of domains fl where separation of variables is possible (explicitly, those which are 
a product of intervals in a coordinate system in which A is separable [2"U]1. 

Many applications demand high eigenfrequency kj (i.e. high mode number j), 
which creates a challenging numerical problem. For instance, the design of high- 
power micro-laser resonators [SS] requires j > 10 3 (i.e. tens of wavelengths across 
the domain). Knowledge of eigenfunctions informs high-frequency wave scattering 
from resonant structures such as jet engine inlets |38j . Interest has also surged 
recently in quantum chaos [621 145] and spectral geometry |28j . where numerical 
studies have played a key role, such as in the discovery of 'scars' of periodic ray 
orbits in chaotic eigenfunctions |30j , and the study of eigenfunction equidistribution 
rates [9]. This can involve computing thousands of modes at up to j ~ 10 6 , 
i.e. hundreds of wavelengths across the domain [60l[9]. The above motivates the 
creation of efficient high frequency numerical methods with controlled errors. 

Existing numerical methods for ((J)-© generally fall into two classes: 

A. Direct discretization of 57 (via finite differences or finite elements [4]), which 
has the advantage that eigenvalues kj are approximated by the spectrum 
of a linear (sparse, often generalized) matrix eigenvalue problem. How- 
ever, since several degrees of freedom per wavelength in each dimension 
are needed, the number of unknowns N grows at least like k d . In fact, to 
achieve bounded accuracy as k — > oo an increasing number of unknowns per 
wavelength are required; this is the so-called 'pollution effect' [5]. Iterative 
methods are needed for such huge eigenvalue problems. We believe the fur- 
thest this has been pushed in d = 2 is j ~ 3 x 10 3 (around 30 wavelengths 
across the domain), by Heuveline and others [3H [SSI [H] . However, here 
specialized multigrid and removal of spurious eigenvalues are needed, and 
relative errors in kj are as high as 10 -3 . 

B. Boundary-based methods, which make use of a basis of analytic solutions 
to the Helmholtz equation (fT|), hence only require discretization of dtt via 
a much smaller N = 0(k d ^ 1 ) unknowns. The main disadvantage is that, 
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since the fc-dependence of the basis is nonlinear, eigenfrequencies kj are now 
given by a (dense) nonlinear eigenvalue problem. This generally requires 
repeated iterative minimization of some error measure along the k axis, 
which is cumbersome and prone to the omission of eigenfrequencies [6l I57j . 
The error measure is often a minimum singular value (e.g. see App. IB)), 
hence 0(N 3 ) effort is required per eigenfrequency found. 

This class includes the method of particular solutions (MPS) [H] (also 
known as collocation, Trefftz, non-polynomial FEM, or ultra-weak vari- 
ational formulation [THl E2]) which uses plane- wave [3D], regular Bessel 
[17l [52] , or corner-adapted Fourier-Bessel solutions [23j [XSJ [14] ; the method 
of fundamental solutions [33] which uses point sources placed outside of f2; 
and boundary integral equation (BIE, also known as boundary element) 
methods which make use of potential theory on <9i7 [19]. Such methods 
often have spectral (i.e. super-algebraic) error convergence, although most 
BIE implementations remain low-order [35] |6] [24] [59]. They can easily 
reach j = 10 4 , with relative errors as small as 10 -14 [12], and variants have 
reached j > 10 6 [60j[57]. 
Can one combine the advantages of classes A and B, i.e. is there a boundary-based 
method that does not require a root search for each eigenfrequency? This was an- 
swered, in the case of star-shaped domains, by Vergini-Saraceno [5]] who proposed 
a 'scaling method' — reviewed in section [7] — which may be viewed as an acceler- 
ation technique for the MPS. Here a single dense matrix eigenvalue problem, i.e. 
effort 0(N 3 ), approximates all eigenfrequencies (and their eigenfunctions) lying in 
an interval of the k axis of length e = O(l). Since by Weyl's law [26] Ch. 11] 
one expects 0(k d ^ 1 ) such eigenfrequencies, this is also the speed-up factor of the 
method, assuming errors are acceptable. The absolute eigenfrequency error is em- 
pirically 0(e 3 ) [3 [9], although this has scarcely been studied. The scaling method 
has allowed large-scale studies of quantum chaos to be performed in d = 2 [60]l9l[TT] 
and d — 3 [48] at speeds around 10 3 times faster than any other known method. 

This key idea of linearizing the nonlinear eigenvalue ■problem in class B has been 
noticed by couple of other researchers. Kirkup-Amini [35] used the linear formu- 
lation of a polynomial eigenvalue problem to approximate the nonlinear eigenvalue 
problem, for low k only. In terms of BIE, Tureci-Schwefel [57] have used the em- 
pirical observation that as a function of k, eigenvalues of the double layer operator 
(see ([30]) below) rotate in the complex plane at roughly constant speed. Veble et al 
convert the BIE to a generalized eigenvalue problem to similar effect [SS]. Heuris- 
tically, these last two methods have the same O(N) acceleration as the scaling 
method. However, the error analysis of the scaling method or such variants is very 
primitive, and certainly no rigorous results exist. 

Here we remedy this by presenting, and analysing in depth, a new class B lin- 
earization method for the eigenproblem (HJ-© in smooth star-shaped domains, 
close in spirit to a BIE method. It is based upon the fc-dependence of the spectrum 
of an intcrioiQ Neumann-to-Dirichlet (NtD) operator for the Helmholtz equation at 
wavenumber k, as presented in section [2] The key idea is that an eigenvalue of the 
NtD reaches zero whenever k reaches an eigenfrequency kj , and thus by computing 



In contrast, it is the exterior NtD or DtN map that plays a common role in applying radiation 
conditions in wave scattering. The interior NtD has been used in analysis of inverse problems, 
|43| and to bound eigenvalues 1251 . 
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all small eigenvalues of NtD one may predict all nearby eigenfrequencies kj. The 
basic algorithm is presented and tested in section [3] 

We devote a large part of this work to the analysis of the spectrum and eigen- 
functions of the NtD, in particular their flow with k, in the k — > oo limit. This 
enables us to analyze the basic method (in section [4|, then propose (section [5]) and 
analyze (section[6]) higher-order accurate variants. The main tools we need are: an- 
alytic perturbation theory (in section [2] and App.|A|, microlocal analysis (App.JCj, 
and a generalization of a recent 'spectral window quasi-orthogonality' result of the 
authors [T2] (App.|D|). 

Here we summarize our main theoretical results: 

• When correctly weighted (as in [5T]) by the function {x ■ n) -1 on dfl, the 
spectrum of the NtD varies approximately linearly with k with slope known 
a priori (Theorem l4.1[) . This will imply that the basic linearization method 
has the eigenfrequency error estimate \kj — kj\ < C(e 2 jk + e 3 ), where kj is 
the approximate eigenfrequency. It is crucial that here C is independent of 
k. Since we establish a one-to-one correspondence between eigenfrequencies 
slightly larger than k and slightly negative NtD eigenvalues, this proves that 
our method has neither spurious nor missing eigenfrequencies. 

• We propose a higher-order accurate formula for prediction of eigenfrequen- 
cies ([6"5"]1 . using an identity for Hclmholtz solutions (jT7|) due to the first 
author 9 . We show that the dominant error term is 0(e 5 ), which im- 
proves upon the 0(e 3 ) of existing scaling methods [TT | [59 ] , 

• We propose a higher-order accurate formula for (boundary data of) eigen- 
functions (|77| . This requires formulae for the 1st and 2nd derivative with 
respect to k of the NtD cigenfunctions k is an eigenfrequency (Prop. [52]). 
We will show a dominant error L 2 -norm of 0(e 3 ). 

• In d = 2, and making a spectral non-concentration assumption (see As- 
sumption [6T]), we prove rigorously that these higher-order methods achieve 
a dominant eigenfrequency error of 0(e ) (Prop. [6~5]) and eigenfunction er- 
ror of 0(ke 3 ) (Prop. 16. 4j) . The latter has the same dependence on e as 
existing scaling methods. 

• In addition, we believe that Lemma 14.41 and the results of Appendix [A] 
are useful contributions to the theory of elliptic boundary- value problems, 
independent of any numerical considerations. 

On the implementation side, we show in section 13.31 that the spectrum of the 
above NtD operator may be approximated with an error uniformly close to machine 
precision, hence that the above error bounds hold in practice. This requires a new 
method based upon potential theory, the Cayley transform, and the quadratures of 
Kress [36]. For <90 an analytic curve, we demonstrate exponential convergence. This 
improves upon the low-order quadratures of all previous scaling variants [61] [57] [59] 
and almost all BIE methods for eigenvalues in the literature. 

We compare the performance of our method in d = 2 against a standard BIE 
root-search (described in App. [Bj, which also serves to give us reference sets of 
kj and <pj against which to measure errors. We test two domains, one with no 
symmetry, and, in section ["5. 31 one with a symmetry that causes an abundance of 
degeneracies. The latter is evidence that Assumption 16 . 1 1 may be violated with no 
impact on performance. We find that the O(N) speed-up translates in practice to 
a factor of 10 3 faster solution at high frequencies. 
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Figure 1. (a) Flow of the eigenvalues /3(fc) of the weighted 
Neumann-to-Dirichlet map vs wavenumber fc, for the domain 
f2 given in polar coordinates by r(9) — 1 + 0.3 cos[3(# + 0.2 sin0)]. 
(b) Eigenmode ^93 (density shows absolute value; white is zero) 
with fc 93 = 19.94995891589 • • • (also shown by red dot in a), (c) 
Zoom in of flow, (d) Zoom in to the same region for eigenval- 
ues of A(fc) -1 the unweighted Neumann-to-Dirichlet map; there is 
variation in slopes at the zero-crossings. 



In section [7] we give a new understanding of the original Vergini-Saraceno scal- 
ing method in a mathematical framework, and draw some comparisons with our 
proposed method. Finally, we conclude in section [S] and give some open questions. 
We have made a documented software implementation of the proposed algorithms 
freely available (in the MPSpack toolbox for MATLAB), and intersperse section [3] 
and beyond with code examples showing how to use these routines. 



2. The Neumann-to-Dirichlet map and its spectral flow 

We will use u n := d n u := n ■ Vit to denote the outward normal derivative of a 
function u defined in Q. Let A(fc) be the interior Dirichlet-to-Neumann operator 
for the Helmholtz equation with parameter fc 2 , that is, the operator that sends a 
function g G iJ 1 (9fi) to the function h £ L 2 (dVt) given by h = u n , where u is the 
interior Helmholtz extension satisfying 

(A + k 2 )u = in u\ dn = g. (3) 

This is well defined for every fc € C except when k — kj is a Dirichlet eigenfrequency, 
in which case the function u may not exist (and is nonunique when it does exist). 
It is well-known that for fc ^ kj, A(fc) is self-adjoint, as the following elementary 
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calculation involving Green's second identity shows. Suppose that g,h £ 
and that u, v are their interior Hclmholtz extensions, respectively, then 



(Ag)h — / gAh = / u n v — / uv n 
an Jan Jan Jan 



= [ [(A + k 2 )u]v-u{A + k 2 )v 
Jn 







This and other properties of A are presented by Friedlander |25j . 

Unless indicated, we work with a weighted inner product on the boundary 8Q, 
denoted by angle brackets (•, •), and induced norm, as follows, 

(g, h) := / gjsjh(s) (x(s) ■ n^))" 1 ds , ||. 9 || 2 := (g,g) . (4) 
Jan 

Note that if f2 is strictly star-shaped about the origin, the weight is bounded and 
positive. It is easy to check that the operator (x-n)oA(k) is self-adjoint with respect 
to this weighted inner product. Let 0(fc) denote the inverse of this operator, that 
is, Q(k) := A(k)~ 1 o (x-n)^ 1 , then O is also self-adjoint with respect to the weighted 
inner product. By definition, if u is any interior Hclmholtz solution, then 

0(fe) (x ■ n)u n = u\ dn . (5) 

In this paper, we will analyze the flow of eigenvalues and eigenspaces of Q(k) as k 
varies along the real axis, that is, nontrivial solutions / G L 2 (dfl) to 

0(fc)/(fc) = (3(k)f(k) ■ (6) 

Taking u = <f)j, a Dirichlet eigenmode, in ([5]), we see that 0(fc) has a zero eigenvalue 
at each Dirichlet eigenfrequency k = kj, with eigenfunction / = (x ■ n)d n <f>f, this is 
why the Neumann-to-Dirichlet map is of interest computationally. Considering the 
case of u a Neumann Laplace eigenmode of the domain shows that 0(fc), and hence 
its spectrum, has a pole at each Neumann eigenfrequency k. Fig. [T] illustrates the 
zeros in /3(fc) occurring at each of the lowest 93 Dirichlet eigenvalues of a domain 
(the poles are also hinted at for larger negative /3). Also visible is the accumulation^ 
of eigenvalues at + that occurs for all real k, a result of A being a pseudodiffcrcntial 
operator of order +1 [25], hence a compact operator (of order —1). 

We wish to flow along an interval of the real k-axis that will likely contain several 
Neumann eigenfrequencies, and need to guarantee that all of the eigenprojections 
and eigenvalues of vary smoothly except possibly for a finite number associated 
with a pole if k is a Neumann eigenfrequency. To do that, we consider the Cayley 
transform of 0, 

R(k) = (0(fc)-i)(0(fc) + i) _1 . (7) 

In Appendix[Xl Corollary |A.2[ we show that R(k) is analytic in some neighbourhood 
of the positive real axis. As R(k) is unitary for real k, its spectrum lies on the unit 
circle, and is discrete except at — 1 because is compact (its spectrum accumulates 
only at 0). Thus we deduce from Kato [34, Ch. VII, sec. 3] that the eigenprojections 
and eigenvalues of R(k) vary analytically away from eigenvalue —1. Translated back 
to this means that the eigenprojections and eigenvalues of 0(/c) vary analytically 
in k on any finite k- interval away from eigenvalue 0, apart from a finite number 
which have a pole at one of the Neumann eigenfrequencies in this interval. 



2 The small gap visible above is due to the numerical approximation of the operator. 
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Definition 2.1. Let f3 = f3(k) be a finite eigenvalue of 0(fc) with boundary- 
normalized eigenfunction / = /(fc), ||/|| = 1. The extended eigenfunction is then 
the unique solution u to the interior boundary-value problem 

(A + k 2 )u = in fl , (8) 
(x ■ n)u n = f on dil , (9) 
u = Pf on dfl . (10) 

Note that we have both Neumann and Dirichlet conditions on it; the latter is 
needed for uniqueness when k is a Neumann eigenfrequency. Their consistency at 
all k is ensured by ([5]). We may view it as a solution to a Stekloff eigenvalue problem 
with Robin condition 

(3(x ■ n)u n — u . (11) 

Note that the extended eigenfunction u is not normalized in L 2 (tt). 

The rate of change with k of each isolated eigenvalue is then given by the following 
variant of a result of Friedlander [25] Prop. 2.5]. For convenience we give the proof. 

Lemma 2.2. Let /3(k) be a analytic eigenvalue branch of 0(fc) with normalized 
eigenfunction f = f(k). \\f\\ = 1, and let u be its extended eigenfunction. Then, 
using the notation (3 :— d(3(k)/dk, it holds that 

(3 = 2k [ \u\ 2 . (12) 
Jn 



Proof. From ([6]) follows the usual Hellman-Feynman formula, 

$ = (©/, /) = (©/ + e/, /) + (e/, /) = (Of, f) (13) 

where the last step comes from the normalization of /, which implies (/, /) = 0. 
Let k = fco be the frequency in the statement of the Lemma, and restrict for now to 
the case that this is not a Neumann eigenfrequency, in which case there is a unique 
solution u to the pair ([5} and © given boundary data /. Holding this boundary 
data fixed at f — /(fco), let v(k) be the solution to the boundary value problem 

(A + fc 2 Mfc) =0, (x ■ n)v n (k) = f(k ). (14) 

(Note v is not the same as the extended eigenfunction u except at k = fco-) Let v 
be the fc-derivative of this solution at fc = fco- Then, by the definition ([SJ), 

©/ = v\an at fc = fc . (15) 

Also, by differentiating the defining conditions (|14p we get a boundary value prob- 
lem for v, 

(A + k 2 )v = -2kv in fi, v n = on dn. (16) 
Combining this with (fl5j) and f| 13[) in Green's 2nd identity gives 



/3(fc ) = (6(fc )/, /) = / (x ■ n)- l vf = / vv n = / (vv n - 
Jan Jan Jan 

1J{A + k 2 )v - [(A + k 2 )ij]v = 2fc / \v\ 2 =2k Q [ \u\ 2 



(17) 



This completes the proof when fco is not a Neumann eigenfrequency. When fco is 
a Neumann eigenfrequency, /(fc) and (3(k) are still analytic in a neighbourhood of 
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ko (as discussed above), so one may take a sequence with ko as the limit and prove 
the same formula. □ 
This fact that this lemma can be applied in the limit j3 f is justified at the end 
of App. [A] Notice that we always have $ > 0, illustrated by the positive slopes in 
Fig.ffl 

We now can explain the reason for choosing the particular weight in the inner 
product (QJ. Let /3(h) be an analytic eigenvalue branch of Q(k) which has _9(kj) = 
for some j, that is, the branch corresponding to the jth eigenfrequencyo Then at 
k = kj, the extended eigenfunction u is a Dirichlet eigenfunction. For Dirichlet 
eigenfunctions, we have Rellich's identity [50] (a special case of (J47J) ) , 

2k] [ \u\ 2 = [ (x ■ n)\u n \ 2 = </,/)= 1 . (18) 
Jn J on 



Inserting this into Lemma 12.21 gives a formula for the slopes at zero eigenvalue, 

(3 = P(kj) = ^. (19) 

Kj 

Remark 2.3. (fT9)) shows that, for the special boundary weight function (x ■ n) -1 , 
the eigenvalues of Q(k) cross zero at a uniform, predictable positive speed that is 
independent of the details of the distribution of the eigenmode <pj . This predictable 
behavior is not known to occur for any other weight function: for example, the 
contrast between this special weight and the unweighted case (where speeds vary 
unpredictably with j) is shown in Fig. [T] (c) and (d). 

3. Basic numerical algorithm 

We first present a simple fast algorithm to approximate the eigenfrequencies 
and eigenfunctions of the domain Q using spectral data of 0; in section [5] we will 
improve it to have higher-order accuracy. 

3.1. Reconstructing eigenfrequencies. Since each Dirichlet eigenfrequency kj 
is associated with an analytic eigenvalue branch of the spectrum of Q(k), we may 
use this spectral flow of Q(k) to locate approximately the kj. Fig.[T](c) illustrates 
that the gradients f3{k) are approximately constant on each branch for k near kj] 
in section |4] we will prove that the range of kj — k for which this usefully holds is 
a constant independent of kj. Thus, choosing a frequency fc* and computing the 
spectrum of 0(/c*), then for each of its small negative eigenvalues /?* = one 
may extrapolate linearly to the corresponding Dirichlet eigenvalue by the 

k 

"linear estimator" : k = — . (20) 

1 + & 

This follows simply from (|19[) and by making the linear approximation f3 t ~ 
f3(kj)(k* — kj). We keep only those k values lying in the interval or 'window' 
[fc*, fc» + e], where e is an O(l) constant. Since, by Weyl's law [25J Ch. 11] asymptot- 
ically 0(k d ^ 1 ) eigenfrequencies lie in such an interval, this is also the order by which 
the method is faster than the standard iterative search for each eigenfrequency. By 
repeating the above with adjacent intervals one may find approximations to all 
eigenfrequencies lying in any desired subset of the frequency axis. 



^This existence of this branch is guaranteed by Proposition I A . 51 
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Figure 2. Errors with basic method, for the domain of Fig. [TJb). 
(a) Error of predicted eigenfrequency kj using linear formula (|20|) . 
vs prediction distance ej, for all frequencies kj £ [90,100]. Lines 
show 0.27e 3 and 0.25e 2 /k. (b) Errors of predicted boundary func- 
tions / in the weighted L 2 norm (QJ. Line shows 0.25e. 

In section 13.31 we present the spectrally-accurate method we use (in d = 2) to 
compute numerically the spectrum of 0(fc»). This algorithm has been built into the 
MPSpack toolbox toolbox for MATL AB [54] , so that the set of approximate eigenfrc- 
quencies kj lying in [90, 100] may be computed, for example, for the nonsymmetric, 
smooth (in fact analytic) domain fi C K 2 shown in Fig. Q] (b), as follows: 
s = segment . smoothnonsym (720, 0.3, 0.2, 3); 7» create a closed curve 
d = domain(s, 1); % create an interior domain 

s.setbc(-l, 'D'); °/ Dirichlet BCs on inside 

p = evp(d); % create eigenvalue problem 

o.khat = '1'; o.eps = 0.1; p . solvespectrum( [90 100], 'ntd' , o) ; 

Here N = 720 sets the number of boundary quadrature points to about 6 per 
wavelength on the boundary, typically sufficient for approximating at close to 
double-precision accuracy. The options structure o chooses the linear method (|20|) 
and sets e = 0.1. The object p now contains p.kj, being the list of 492 approxi- 
mate eigenfrequencies kj found (these are in fact numbers j — [2064, 2555] for the 
domain). All were found to be simple, as expected generically since Q has no sym- 
metry. The majority of them have absolute errors less than 10~ 4 . The CPU time 
for the above example was 13 min, ie 1.6 s per computed eigenfrequency^ 

The size of the absolute eigenfrequency errors are shown in Fig. [2] (a), versus 

Cj . — kj k* , (21) 

the frequency 'distance' over which the linearization occurred. Errors are 0(e 2 /k) 
at small distances but O(e^) at large distances: these two terms are shown by 
straight lines in Fig. [2] (a) . We are able to prove a bound involving these two terms 
in Corollary 14.21 which states that the implied constants are independent of k. The 

^Runtimes are reported for a 2005-era workstation (two single-core Opteron 2GHz 250 CPUs) 
with 8 GB of RAM, running linux, MATLAB 2008a, and MPSpack version 1.2. 
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transition point (intersection of the straight lines) occurs at e = 0(l/k). Since 
generically only a fraction 0(1/ k) of the eigenfrequencies in the window lie below 
this e distance, the method is asymptotically 3rd-order accurate in the interval 
width e. 

These errors reported above were found by comparison against an accurate set of 
eigenfrequencies kj found independently by a standard method from the literature 
described in App.JEJ This reference method requires 53 s per eigenfrequency found, 
thus our method is a factor 33 times faster. Assuming constant absolute eigenfre- 
quency error is acceptable, then this speed-up factor grows (in d = 2) in proportion 
to O(N) — 0(k): the reference method takes 0(N 3 ) effort per eigenfrequency found 
whereas our proposed method takes only 0(N 2 ) effort. 

3.2. Reconstructing eigenfunctions from boundary data. We assume for 
now that for each eigenvalue /3* of 0(fc») we can generate an accurate approxima- 
tion to its corresponding boundary eigenfunction /* (e.g. as in section [33)1 . Ap- 
proximations 4> to Dirichlet eigenfunctions cf>j can then be evaluated using potential 
theory, as follows. 

At wavenumber k, the free space Green's function for the Helmholtz equation, 
Go(k; x, y), is defined as the unique radiating solution to — (A + k 2 )Go = S in R d , 
where 5 is the Dirac delta distribution. Specifically, we have, 

. / , \ rf/2-1 

G (fc; g ,y) = l ^-j— — J H^iklx-y]), x,yeR d , (22) 

where H„ is the outgoing Hankel function of order v 46, Ch. 10]. The standard 
single- and double-layer potentials [18] are then defined for x € f2 by 

(S(k)a)(x) = f G Q (k;x,y)a(y)ds y , (23) 
Jan 

(V(k)r)(x) = f dG f> X > y K (y)ds v , (24) 
Jon on y 

where the derivative is with respect to the y variable in the outward surface normal 
direction at y. Then any solution m to (A + k 2 )u — in fl with smooth boundary 
may be written via Green's representation theorem [18] , 

u = S(k)u n — T>(k)u\gii in Q . (25) 

Suppose an exact eigenfrequency kj were known, and also fj = f(kj) the cor- 
responding exact eigenfunction of Q(kj). We could then use ({25]) to compute the 
extended eigenfunction Uj, since its Dirichlet data vanishes, and its Neumann data 
is given by ((9|) . According to (|18|) we also need a prefactor <f>j — \/2kj Uj to re- 
cover unit L 2 (Q) norm. Thus a Dirichlet eigenfunction <pj is represented exactly 
throughout £1 by 

0j = s/2kj S(kj)[(x ■ ny 1 ^] . (26) 

However, we do not have access to fj\ we only have /*, the corresponding eigen- 
function of 0(/c*) for fc* near kj. Similarly, kj is only known approximately (e.g. 
as in the previous section). Given approximations k kj and / ~ fj t we propose 
to reconstruct an approximate eigenfunction via 

4> = V^kSfyftx-n)- 1 /] . (27) 
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\\fj-fj\\ 

Figure 3. Scatter plot of estimated eigenfunction errors in the 
L 2 (f2) norm versus corresponding boundary function / errors in the 
weighted L 2 (dfl) norm, for the domain of Fig. Q] and kj € [90, 100] 
with the basic method. L 2 (f2) norms are estimated on a Cartesian 
grid of 327 interior points, giving a statistical error of order ±10%. 

For now we will present a method that is only first-order in e: we use the 

"trivial / estimator" : f = f* ■ (28) 

(We present higher-order methods in section [5]) Figure EJb) shows the resulting 
||/ — /ill errors in the weighted L 2 (dfl) norm, computed relative to a highly-accurate 
set of boundary functions /» found by the method of App.lBl This behavior is clearly 
first order. 

Remark 3.1. To prove a rigorous estimate on ||/ — fj\\ one would need to control 
/ over the interval we have by ([71)1 and Lemma [4.41 that ||/|| = O(l) at 

k = kj, but cannot exclude the possibility that "avoided crossings" in the spectral 
flow cause / to be much larger at other k values. Based on empirical observations, 
the latter possibility seems very rare. 

How do the errors in / propagate to errors in eigenfunctions cj>? To test this, we 
insert / = /*, and k from (|2"0"|) . into the reconstruction formula (|2"T1) . and estimate 
numerically the L 2 (Q) errors against an accurate set of reference eigenfunctions 
4>j . In the resulting Fig. |3] the data clusters close to a straight line of unit slope 
(scatter from this line being part due to our estimation of L 2 (Q) errors using a 
relatively small number of interior points). Hence the domain error norm of <fi is 
empirically controlled by the boundary error norm of /. This is to be expected 
because, although (|26p and (|2"T)) use different wavenumbers, the error in fc is of 
higher order than that of /, and error induced by the fc-dependence of S(k) is 
expected to be negligible. 

Remark 3.2. Supported by the above evidence, we henceforth discuss eigenfunction 
errors only in terms of boundary functions /, postponing analysis of \\<p — <fij ||l 2 (s"2) 
to future work. A rigorous proof that boundary error controls domain error would 
demand bounds on the fc-dependence of the operator S(k) : L 2 (dQ) — > L 2 {VL). 
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Figure 4. Convergence of numerical scheme for eigenvalues j3 of 
0(fc) at three different k values with distances from a Neumann 
eigenfrequency as follows: (a) 10~ 2 , (b) 1CT 10 , (c) zero. The do- 
main is as in Fig. [TJ The (3 tested was the largest negative eigen- 
value (i.e. closest to zero), roughly —5 x 1CP 3 in each case. The 
proposed Cayley scheme ([55)1 and ([57)1 is compared against the 
direct discretization of ([551) . 

Accurate numerical study of L 2 (f2) errors is also difficult, since i) the eigenmodes 
are highly oscillatory, demanding 0(k 2 ) evaluation points (in the above example 
around 10 5 would be needed), and ii) accurate evaluation of a layer potential such 
as ([27| near Oft is difficult and a topic of current research [31] . 

In terms of computational effort, extracting all boundary eigenfunctions /* at 
each fc* is best done by complete diagonalization of a matrix (given below by (|44|) ) 
at each A;*; this increases the CPU time per mode from the 1.6 s of the previous 
section (when only matrix eigenvalues were needed) to around 2.3 s per mode. 
However, the reference method of App. [Blalso requires longer to extract modes (an 
additional 14 s per mode). The net effect is that the proposed NtD method is still 
30 times faster than the reference method. 

Remark 3.3. In [12] we proved error bounds on k and on the L 2 (fl) error of (ft 
in terms of \\(p\\ l 2 (9H)- The latter could be evaluated using (f27| and a singular 
quadrature scheme as in Section 13.31 This would remove any need to compare 
against reference eigenpairs. However, we avoided this approach since the errors in 
/ would dwarf the higher-order errors in k that we wish to study. 

3.3. Numerical computation of spectrum of 0. To implement the above al- 
gorithm, at any given frequency k we need to compute numerical approximations 
to an O(l) fraction of the eigenpairs of the weighted NtD operator 0(fc). Here 
we present, and test, a robust integral equation method based upon the Cayley 
transform. We need some standard results from potential theory [18] , Let S(k) 
and D(k) be the single- and double-layer boundary integral operators formed by 
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restricting ([23]) and (|24|) respectively to the boundary, 

(S(k)a)(x) = [ G (k;x,y)a(y)ds y , x £ dQ (29) 

(D{k)r){x) = f dG °^ Xiy K {y)ds y . x £ dQ (30) 
Jdn (J n v 

Then, taking x £ fl to the boundary in the representation formula (l25j) , and apply- 
ing the jump relation for the double layer potential, 

{Vr)\ aQ = {D-\)T , (31) 

gives, for any Helmholtz solution (A + k 2 )u — in O, the boundary data relation 

{\ + D)u\ aa = Su n . (32) 

We also record for later use the jump relation for the single layer potential, 

(Sa) n = (D« + i)a . (33) 

We now generalize the Cayley transform ([7]) slightly, defining 

R ri {k) := ( V e(k) - i) (r,e(k) + i)' 1 , (34) 

where 77 £ K \ {0} is a scale parameter with units of inverse length (i.e. of k). 
Therefore, if R v (k)f — g, we have 

( V e + i)g=(r)G-i)f, 

which we rearrange to 

®[in{g-f)] = 9 + f- 

That is, there exists a function u on Q, such that 

(A + /s 2 )u = 0, u n = i-q{x ■ n) _1 (.g - /), u\ dn = g + f. (35) 

Inserting this boundary data into (1321) implies 

i ^[(x-n)- 1 ( ff -/)] = (i+^)(.g + /) 

which can be rearranged, recalling that R v f = g for all /, to show, 

R V = {K^)- 1 K+, K± = ±{\ + D) + ir)S o {x ■ n)- 1 . (36) 

The scheme is now to choose the scale parameter (we prefer r\ — k, similar to [36]), 
and to approximate the spectrum and eigenfunctions of i? r) , using known efficient 
Nystrom discretizations for the operators S and D, as described below. We then 
convert back to eigenpairs of as follows: the eigenvalues (3 of come from the 
eigenvalues A of R n simply by inverting the formula (l34l) . that is, 

P = ~, (37) 
7] 1 — A 

and the eigenfunctions of are the same as those of R v . 

Remark 3.4. The advantage of discretizing (|36[) then transforming eigenvalues via 
(|37p , over discretizing the usual direct representation of the weighted NtD map 

6= {± + D)- 1 S o (x ■ n)- 1 (38) 

which follows from (|32[) . is that Rr){k) is unitary and thus its eigenvalues remain 
of size O(l). By contrast, the eigenvalues of 0(fc) have a large dynamic range, 
and a finite number of eigenvalues diverge to infinity whenever k is a Neumann 
eigenfrequency of the domain, causing inevitable large round-off error in the desired 
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(small) eigenvalues. We demonstrate this contrast numerically in Fig. 2) in the 
'direct' method this round-off error limits accuracy in /3 to 8 digits for k near a 
Neumann eigenfrequency (and fails to converge at a Neumann eigenfrequency), 
whereas the 'Cayley' method achieves 14-digit accuracy uniformly in k. (Note that 
we expect some mild loss of accuracy as k increases, due to the condition number 
of the (K-)^ 1 factor, but in the range explored in this paper, 1 < k < 10 3 , this 
was negligible.) Thus to discretize (|38|) is not robust, whereas our proposed scheme 
is robust. 

We summarize briefly our preferred Nystrom discretization for Helmholtz layer 
potential operators on analytic curves in d = 2, following Kress [35]. Let z : 
[0, 2ir) — > M 2 be a 27r-periodic parametrization of dil, and let k(x, y) be the kernel 
of either S or D. Changing variable to s,t G [0,2-71") we get kernel K(s,t) :— 
k(z(s), z(t))\z'(t)\ where z' = dz/dt. Note that S has a logarithmic singularity on 
its diagonal, whereas D has a continuous kernel but is non-analytic on the diagonal; 
in both cases the following splitting allows spectral accuracy to be achieved. We 
choose quadrature nodes tj — 2nj/N, j = 1, . . . , N, and split the kernel into the 
form 

K(s,t) = log ^4 sin 2 tZl^j K 1 {s,t) + K 2 (s,t) (39) 

with K\ and K% both 27r-periodic and analytic. The matrix representation of K 2 
comes from the periodic trapezoid rule (weights being constant at 2ir/N), whereas 
the representation of K\ involves a product quadrature appropriate for the peri- 
odized log singularity. Together these give a matrix K with elements 

2tt r (JV) 



K 



13 ~ N 



R\^(0) K^t,) + K 2 (U,t 3 )\ , i,j = l,...,N, (40) 

where the Martensen-Kussmaul quadrature weights (deriving from the Fourier se- 
ries for the log factor, see jHTJ Lemma 8.21]) are defined by 

W /2-l 

R^is) = - V -cosm(s-t 7 ) -cos— (s-tA . (41) 

1 w ^ m N 2 v ; 

m— 1 

Abusing notation slightly by letting K be an operator with kernel K(s,t), it is 
standard to approximate operator eigenvalue problems of the type 

K4> = \<t> (42) 

by the TV-dimensional matrix eigenvalue problem 

= xWjW . (43) 

If K were compact and normal, it is known that the spectrum and eigenspaces of 
(l43l) converge to those of (|42|) as N — > oo [2] , at a rate given by the error of the quad- 
rature scheme applied to vectors in the eigenspace (for the spectrum see [3] — here 
normality ensures that the index of each eigenvalue is one — and for the eigenspaces 
see [271 Thm. 1]). This analysis relies on the framework of collectively compact 
operators pQ |37l Ch. 10]. The above product quadrature scheme is within this 
framework and is spectrally accurate for analytic functions, i.e. errors are bounded 
by ce~^ N for some 7 > [36l [37]- 

However, our goal is to approximate the spectrum and eigenspaces of the operator 
R n : this is not covered by the above-mentioned analytic results, for two reasons. 
Firstly R v is not compact (although R v + 1 is), and secondly the application of R v 
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in (|36|) requires an operator product and inverse. We will not attempt a rigorous 
error analysis here, rather merely describe our scheme and show its efficacy. We 
approximate the spectrum of by that of the matrix 



built from the matrices K± which approximate the operator factors K± appearing 
in (|36l) . according to the above Nystrom scheme (I4U1) . Dense linear algebra is 
used both for the matrix inverse KZ 1 in (|36|) . and the full diagonalization of 
(MATLAB's inv and eig respectively). The computational effort is 0(N 3 ). The 
eigenvectors of R r) then give approximations to the eigenvectors of R^, hence of O, 
at the quadrature nodes. In MPSpack the above algorithm is available via 
[beta.V] = p.NtDspectrum(k) ; 

which returns approximate eigenvalues of 0(fc) in beta, and corresponding eigen- 
function values at the quadrature points in the columns of V. Returning to Fig. |4]we 
observe exponential convergence of this 'Cayley' scheme, with saturation at relative 
error 10~ 14 uniformly in fc. 

4. Error analysis of the linear eigenfrequency estimator 

The main result of this section is an estimate on the accuracy of the eigenfre- 
quencies as reconstructed by the basic formula ([20)1 . In section I5TT1 we will describe 
an improved method for which we can prove better error estimates, but those better 
estimates are conditional on absence of spectral concentration (Assumption 16.11) ; 
here, the result is unconditional. 

The key result is the following bound on the deviation from linearity of the 
weighted NtD eigenvalue flow. 

Theorem 4.1. There are constants e, K > dependent only on fi such that the 
following holds. Let k p > K be a Dirichlet eigenfrequency, and fi{k) be the corre- 
sponding eigenvalue branch ofQ(k), i.e. such that /3(k p ) = (the existence of which 
is guaranteed by Proposition \A . 5\ ). Then 



for all k € [k p — e, k p \. The implied constant in the O(-) depends only on fl. 

It is then easy to derive the following error estimate for the basic method. Note 
that we have already numerical evidence (section l3.ll) that the powers of e are sharp. 

Corollary 4.2. There are constants e, C > depending only on Q such that, for 
any sufficiently large fc» 7 and any k p lying in the range [fc*, fc* + e], and /3(k) related 
to k p as in Theorem \J^.l\ we have 



where k — fc*/(l + /?(&*)) is the linear estimator of k p according to (|2Q|) . 

Remark 4.3. Note that the theorem holds for a fixed window width e, independent 
of fc. By Weyl's Law there are 0(k d ~ 1 ) eigenfrequencies lying in such a window; 
all are found within the stated error. As fc grows, the e 2 /fc term becomes negligible 
for almost all reconstructed eigenfrequencies, thus the eigenfrequency error of the 
basic method is effectively 0(e 3 ) with constant independent of fc. 



(44) 




(45) 




(46) 
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Proof of Theorem \4-l\ We use the following identity from [9j Lemma 3.1], which 
allows one to express the right hand side of (|T2|) in terms of boundary data: for 
any Helmholtz solution u at frequency k, we have 



2k 2 \u\ 2 = {x-n)(\u n \ 2 + k 2 \u\ - |V tan u| 2 ) + u n Wu + (Wu)T^. (47) 
Jn Jan 

Here, V tan is the tangential gradient on dfl, and W is the tangential part of the 
vector field x ■ V which generates dilations. Explicitly, W — x tan ■ V tan where 
a^tan = x — (x ■ n)n. For example, in d = 2 we have V tan =3t and W = (x • t)dt, 
where t is the unit tangent vector. Putting (|4T|) together with (fl"2"|) . taking u to be 
the extended eigenfunction, we obtain 

^ (fc) "Kij ( x ' n KKI 2 + A:2 M 2 - |V tan w| 2 ) +2Re(«»Wti)V (48) 
The principal term on the right hand side of (|48l) is fusing [TBI . 

\ j { X . n )\u n \^\\\ff = \. 
k Jon k k 

The other terms are small when fj is small, and we try to estimate them in terms 
of f3. Two of the terms are not hard to estimate: we have using the boundary 
condition ([TTj) . 

f (x-n)k 2 \u\ 2 = 0(p 2 k 2 ), 
Jan 

while (using the divergence theorem on <9f2 in the third step below), 

2Re{{Wu)W l ) = [ ix-ny 1 ^ 1 -2Re((Wu)u) 
an Jan 



(x-nj^W^uY) 
an 



\u\ 2 {W + div W)(x-n)- 1 

ian 

= -0 I ((x-n)W{{x-n)- 1 )+ div W)(x-n)\u n \ 2 

Jan 
= 0((3). 

The scalar function div W may also be written V tan ■x taxl . To deal with the |V tall u| 2 
term in (|48|) , we prove the following in Appendix [Cj 

Lemma 4.4. There are constants c, K > depending only on O, such that when- 
ever k > K and u solves (A + k 2 )u — in fl, with 

u = j3(x ■ n)u n on dQ (49) 

for some Robin constant (3 e [— c, 0], then 

||VtW"||z,2(afi) < 2fc||u|| L 2(af2)- (50) 

Remark 4.5. The intuition behind Lemma 14.41 is that /, as a boundary trace of 
a Helmholtz solution at frequency k, should be band-limited to frequencies < k. 
Indeed, the coefficient 2 in ([50]) could be replaced by any factor a strictly larger than 
1, for k > K{a). Also, using the same proof is not hard to show the corresponding 
result for higher derivatives: 

UV&IU^n) < 2k m \\u\\ L 2 (an), k > K m . (51) 
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Using Lemma f4.4[ we can estimate the |Vt an M| 2 term in (|47|) the same way as 
the k 2 u 2 term. So, combining the estimates of terms in (|48|) . we get 

P{k) = I +Or-^) + 0(f3 2 k) whenever - c < f3 < , (52) 



k V k 

with implied constants depending only on Q. 

We now conclude the proof of Theorem 14.11 by establishing (|45|) . This follows 
directly from ()54j) below by integrating in k. Therefore, it remains to prove the 
following result: 

Lemma 4.6. There exists constants e, C\ > depending only on such that, for 
any sufficiently large k p (with (3{k) as in Theorem ^. lty , it holds for all k £ [k p — e, k p ] 
that 

^^<m<^ , and (53) 

Jvp £i iVp 

m-v <^(^ + ^V^) ■ (54) 

Proof. We first prove the left hand side of (|53p . We first use the continuity of 
j3(k) to observe that in some small left neighbourhood (k",k p ] of k p , f3(k) itself is 
arbitrarily close to — in particular, such that f3(k) > — c. Therefore, ([52]) applies, 
and, by requiring j3 sufficiently small we can make the right hand side of (|52|) less 
than 3/ (2k), and hence less than 2/k p (since k/k p can be made as close as we like 
to 1). By integrating this, we find that in this neighbourhood we have the left hand 
inequality in (|53|) . However, the size of the neighbourhood may still depend on k p . 

Now we prove that for some e > 0, the left hand inequality in (I5U1) holds on 
each interval [k p — e,k p ] for all sufficiently large k p . We do this by contradiction. 
Suppose that there is a k £ [k p — e, k p ] such that j3(k) < 2(k — k p ) / k p . Let kl be the 
largest such element of the interval [k p — e, k p ]; notice that k' is strictly less than 
k p using the paragraph above. Then we have 

P(k)>2^^foTk£[k',k p ], P(k') = 2^^. (55) 

K p K p 

For small e relative to k this certainly implies that (3 > — c on the interval [k 1 , k p ], 
so (|52p applies. Using (|5"2")l and the estimate (|5"5")) we find that 

f\ \ ixiivp p 



<i(i + c( 2£ 



1/* \ \ Is* I/ 9 2 

fx \ \ rvp il n 



3_ 

2k 
1.6 

k p 



p (56) 



where we need e sufficiently small in the second last line, and k p sufficiently large 
in the last. Integrating this we find that 



(3(k p ) - /3(k r ) < 



1.6(fc p - k') 
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which contradicts the second part of ([55]). We conclude that no such k' exists, so 
the left hand inequality of ([53| holds on the whole interval [k p — e, k p ] . 

The right hand inequality is proved similarly. In fact, using the left hand in- 
equality, we see for small e that ft > — c on the whole interval [k p — e, k p ], so we can 
use (|52|) on the whole interval, and conclude in a similar way to (|56|) that 

8 

ft( k ) > ~j— , ke[k p -e,k p ] 

to derive the right hand inequality in (|53j) . 

To prove (|54|). we first note that (|53l) inserted into ([52"T) implies that 

m-l\<Ci(^ + ^^), ke[k p -e,k p] . (57) 

This is almost the same as (j5"4")) . but there are factors of k in place of k p . For the 
left hand side, replacing 1/k by l/k p makes an error of (k p — k)/kk p , and this can 
be absorbed on the right hand side (by increasing C\ by 1). Then, by taking k p 
large relative to e, we can replace the occurrences of k on the right hand side by 
k p , at the cost of increasing C\ slightly. We conclude that ([54")) holds. □ 

Remark 4.7. Theorem l4.1l and Corollarv l4.2l inrplv that every k p slightly bigger than 
k* corresponds to a slightly negative eigenvalue ft of 0(£;*), with an almost linear 
relationship between k p and ft. The converse is also true: every slightly negative 
eigenvalue ft of 0(fc*) corresponds to a k p slightly bigger than fc*. To see this we 
note that (f5"l?)) . and Proposition IA.51 imply that the eigenvalue branch starting at 
ft(k*) > —e/k will, for small e, reach zero near k ps fc*/(l + ft(k*)). Thus there is a 
one-to-one correspondence between eigenfrequencies k p slightly bigger than fc* , and 
the slightly negative eigenvalues of 0(fc*). 



5. Higher-order accurate reconstruction methods 

5.1. Higher-order eigenfrequency approximation. In section EO we presented 
a formula (|20[) for eigenfrequencies k p . For its error analysis in section [4] we treated 
all terms on the right hand side of (l48l) . other than the first, as error terms, and esti- 
mated them. However, numerically we have at our disposal not just the eigenvalues 
of 0(fc*), but the eigenfunctions. Observe that the RHS of (|4"5|) can be expressed ex- 
actly in terms of ft and its associated eigenfunction /, using the relation u\gQ = ftf 
from (fTU|). Precisely, we have 

k 2 [ (x-n)\u\ 2 = k 2 ft 2 \\(x-n)f\\ 2 m , (58) 
f (x-n)\V tan u\ 2 = ft 2 \\(x-n)V tlia f\\ 2 dn , (59) 

( 2{x ■ n) Re{u n Wu) = -ft [ (W + div W)((x ■ n)- 1 )]/] 2 =-ft{f, mf) , (60) 
Jon Jen 
where we introduce the scalar boundary function 

m := (x ■ n)W((x ■ n)' 1 ) +divW. (61) 

Using the above, we can rewrite (|4"5|) as 

ft = \ + kft 2 \\{x- n)ff m n)V tao f\\l Q - mf). (62) 
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Figure 5. Errors of the two eigenfrequency prediction schemes for 
four different k ranges, vs ej := kj — k*. Linear scheme (dots), using 
(I2TJ)) . compared against power laws 0.27e 3 and 0.25e 2 /fc (dotted 
lines). Higher-order Riccati scheme (crosses), using (|65[) . compared 
against 0.02e 5 and 7e 3 /fc 2 (solid lines). For the power laws, the 
mean k in the interval is used. The domain is as in Fig. [ljb). 



Of course, for values of k other than fc* we no longer know the exact boundary 
eigenfunction f(k). However, we can get a potentially more accurate estimate of the 
function /3(k) by "freezing" the values of the norms in (|58[) -([60 |) by fixing / = /(&*). 
That is, we define constants 

A:=kl\\(x-n)ff dn -\\(x-n)V tlm f\\l n , B :=-</, mf) , (63) 

where k z is a frozen value of k yet to be specified, and consider the ODE 

P= ±-(l + Al3 2 + B^ . (64) 

After changing independent variable to log A;, this is a constant-coefficient Riccati 
equation that can be solved exactly. Assuming that A > (B/2) 2 which is expected 
(cf. Remark |4.5| note that B = 0(1) as k — > oo), the general solution is 

P(k) = — + -^tan(/zlogfc + a), where fi = y/A- (B/2) 2 , 
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time / mode (sec) 


abs error of kj 


L 2 -error of fj 


k interval 


j 


N 




ref NtD ratio 


max median 


max median 


[30, 40] 


4e2 


300 


176 


8.1 0.72 11 


1.5e-7 1.3e-8 


1.6e-3 1.5e-4 


[90, 100] 


2.6e3 


720 


492 


67 2.3 30 


8e-8 1.2e-9 


3e-3 1.2e-4 


[300, 302] 


2.3e4 


2200 


314 


1200 15 80 


2e-7 3e-~9 


7e-3 3e-4 


[1000, 1000.1] 


2.6e5 


7200 


53 


3e4* 134 250* 


2e-7 6e-9 


l.le-2 6e-4 



Table 1 . Runtimes (per mode found) and errors for proposed NtD 
method for cigcnmodes of the nonsymmetric domain of Fig. [lib). 
The number of modes found in each frequency interval is n m , and 
the approximate mode number is j. In all cases e = 0.1; note that 
error can of course be reduced by reducing e. The Riccati (|65[) 
and quadratic (|77[) estimators were used. In the reference method 
(App. [Bj the absolute kj errors were better than 10 -12 . Asterisk 
(*) indicates estimated values; in fact a faster method (that we 
shall not describe here) was used for the highest reference set. 



where a is an arbitrary constant chosen so that the initial condition /3(fc*) = /3* is 
satisfied. Solving for k p gives the 

"Riccati estimator" : k = fc* cxp — ( tan -1 (- — ) — tan -1 (— )^ . (65) 

Figure [5] shows the observed errors for this Riccati estimator in d = 2 (in our 
code example this is achieved via option o.khat = ' r'); they are 10 3 to 10 5 times 
better than those shown for the linear estimator (|20|) . We in fact compared the 
constant choice k z = fc* against k z = |(1 + (1 + (3*)~~ 1 )k*, the mean of fc* and the 
linear estimator k, and found that the latter choice has slightly smaller errors, hence 
prefer it. We study four frequency ranges, so that the k behavior of the constants 
in the e power laws becomes visible. This provides strong evidence that the error 
of the Riccati scheme is 0(e 3 /k 2 + e 5 ), which is dominated by the second term 
when the window e is chosen to be large enough to collect many eigenfrequencies 
(i.e. > k^ 1 ). Note that the constant in 0(e 5 ) is independent of k, and appears 
quite small, thus absolute kj errors are around 10 -7 for e = 0.1. In Section [6l we 
shall give a theoretical analysis of this method (with the choice k z — fc» ) , under a 
spectral nonconcentration assumption (see Assumption 16. 1[) for 0(k*) at /3*. 

Table [T] summarizes the numerical experiments: note that the speed-up ratio 
relative to the reference method is roughly linear in k, reaching a couple of hundred 
for our largest calculation (around 400 wavelengths across). Thus the speed-up is 
close to the number of wavelengths across the domain, for the errors reported. 

Remark 5.1. We have tested the above Riccati estimator against a more accu- 
rate approximation which considers a linear approximation in k for the quantities 
\\{x ■ n)f\\g n , (|59|) and (|60|) . and solves (|64]l with k varying (this requires a numerical 
ODE solver). We found no significant improvement, hence recommend (|65[) . 
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Figure 6. Boundary error norms of three eigenfunction predic- 
tion schemes for four different k ranges, vs tj :— kj — fc*. Trivial 
scheme (dots), using (|28|) . is compared against 0.25e (dotted line). 
Linearized scheme (circles), using (fT5j) is shown. Quadratic scheme 
(crosses), using (|77|) . is compared against power laws 0.5e 2 fc -1 / 2 
and O.le 3 ^ 1 / 2 . The domain is as in Fig. [TJ 



5.2. Higher-order reconstruction of eigenfunctions. In order to find higher 
order estimators for the Dirichlet eigenfunction (or more precisely, its normal de- 
rivative at the boundary), we first compute the ^-derivative of an eigenfunction 
branch f(k) of 0(fc). For simplicity we assume the eigenspace is simplqj. 
Taking the derivative of ^ gives the formula 

(9 -/?)/=(/?- 6)/. (66) 

Now fix ko £ [fc*, fc* + e] and let v(k) be as in (fl4|) . hence satisfying (|T5j) and (TTBT) . 
Also, at k — k we have 

v(ko)\on=l3(ko)f{k ). (67) 
We make the observation that, due to the commutation formula 

[A,x ■ V] = 2A, 



Note that our rigorous results also make this assumption as it is a consequence of Assump- 
tion ED 
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the function x ■ Vw/fc satisfies 

(A + k 2 ) Q-x ■ Vv^j = -2kv. 

Therefore, combining this with (|16[) . the function q :— v — ^x ■ Vw is Helmholtz 
for every k, and so we have, at k — ko, a relation between the value and normal 
derivative of q on the boundary, 

(0 - /3)(x ■ n)q n = q\ d n - /3(x ■ n)q n . 

Using the second part of (fT6|) this simplifies to the following at k = ko, 



-(6 - 0)(x ■ n)d n (x • V«) = -v\an + hx ■ Vv)\ an - ~(x ■ n)d n (x • Vu) . (68) 

We need to re-express the spatial 2nd-derivatives in terms of the boundary dfl. The 
Laplace-Beltrami operator Agn is related to the Laplacian in K rf by 

A = d nn + (d- l)Hd n + A dn , (69) 

where the scalar function H is the mean curvature of dfl. Thus, for any Helmholtz 
function w, writing the scalar function h := — (d — l)(x ■ n)H , we have, 

((x ■ n)d n ) 2 w\ dn = -{x ■ n) 2 {A an + k 2 )w\an + (I + h)(x ■ n)w n . (70) 

Writing x ■ V at the boundary as (x • n)d n + W, we also compute, for any smooth 
function w, that, 

[(x ■ n)d n ,x ■ V]w| an = g{x ■ n)w n + W'w , (71) 

where g := — (x ■ n)~ l W(x ■ n) is a scalar function, and W' is a tangential derivative 
operator on dtt whose vector field is given by the covariant derivative of n with 
respect to the dilation vector field x • V. Here we extend the normal vector n = 
(ni, . . . , rid) = X) n i e i m to a neighbourhood of dfl so that it is of unit length and 
constant along lines perpendicular to the boundary. Explicitly, 

W':=(x-n) x ^ d s ■ ( 72 ) 

One may check that g + h + 1 = m from (|6Tj) . via the identity h = div W — 1. Thus 
combining (f70| and (f7"Tj) . and inserting (|67| and (x ■ n)v n — /, we get 

(x ■ n)d n {x ■ Vv) = (W + m)f-P(x-n) 2 {A 9n + k 2 )f + f3W'f. 

Acting on this with i(G — /J) then equating with (|68p . replacing v via (fT5]i . and 
again expanding x ■ gives, at k — fco, 

^(9 - 0) ({W + m)f + m'f ~ P(x ■ n) 2 (A d n + k 2 )f) 

= -e/+^(/ + pWf) -^((W + m)f-(3(x-n) 2 (A m + k 2 )f) 

= 08 - 6)/ + (i -/?)/- ^ (m/ + - P{x ■ n) 2 (A dn + fe 2 )/) . 
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Notice that the flWf/k terms canceled in the last step. Combined with (|66|) . and 
observing that the range of — /3 is orthogonal to /, we get 

((W + m)f + (3W'f - f3(x ■ n) 2 (A dn + k 2 )f) 

= (9 - p)f - ^ (rnf + (3W'f - f3(x ■ n) 2 (A 9n + fc 2 )/) ± 

where - 1 indicates projection onto the space orthogonal to /. Now applying (0 — 
/3) _1 (again we consider the generalized inverse, equal to zero on the span of / and 
inverting — /3 on the orthogonal complement), we find 

f=l ((W + m)f + pW'f - P(x ■ n) 2 (A dn + fc 2 )/) 

+ f (6 - (3)- 1 (mf + pW'f - (3(x ■ n) 2 (A dn + k 2 )f) + cf, (73) 

where the constant c is determined by the normalization, i.e. (/, /) = 0. 

From this we can determine the first and second derivatives of f p (k), the eigen- 
function on the branch corresponding to Dirichlet eigenfrequency k p , when k = k p , 
that is, when f3 — 0: 

Proposition 5.2. Let (f(k),/3(k)) be an eigenpair for 0(fc), and let D be the 

differential operator 

D:=W + m-±(mf,f) . 
Then if (3(k) = 0, the first and second derivatives for the eigenfunction f(k) are 

/ = ^((D 2 - D)f - (x ■ n) 2 (A dn + k 2 )f + W'f + Q(k)~ 1 (mffj (74) 

-p(m/,D/)/ + ci/ , 
where C\ is some normalization constant. 

Proof. The first identity follows from ([73)) by setting f3 — 0, and computing that at 
j3 = 0, noting that W + m is adjoint to — W with respect to |4|, 

The second formula follows by taking the /c-derivative of the right hand side of ([73")) 
and then setting f3 to zero, using (|T9| . □ 

This proposition suggests that the following two estimators for f(k p ) should be 
more accurate than the trivial estimator / = f p (k*) considered in Section T3. 2 1 First, 
using just the first derivative formula in (I74[) . we consider, with / = / p (fc*), the 

"linear / estimator": f P -=f + j-Df , where e p := k p — fe* , (75) 

e p being the best available estimate for k p — fc*, e.g. via (165[) . Numerically in d = 2 
we handle the term Wf = (x ■ t)dtf using a N x N spectral differentiation matrix 
[55"! Ch. 3] applied to the discretized /; the FFT could also be used. Referring 
to the data shown by circles in Fig. El for the domain of Fig. QJb), we see that 
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empirically, this estimator is second-order accurate in e, with a constant that is 
independent of k. This improves upon the trivial estimator by one to three extra 
digits of accuracy. 

In principle, one should be able to use the second derivative of / given by (|74|) 
to obtain a third order accurate estimator. Unfortunately, this formula involves 
the operator Q(kp)^ 1 which is not known explicitly; it could be approximated 
numerically at a cost of 0(N 3 ), but this would need to be done afresh at each 
eigenfrequency and thus destroy the 0(N 2 ) complexity per mode. However, if we 
study the size of the terms in the second derivative formula ([71)1 . we see that some of 
them can be expected to be lower order (in k) than others. For example, the terms 
k~ 2 Df and W'f are lower order than k~ 2 D 2 f. Also, as discussed in Remark I D. 3 [ 
subject to a spectral nonconcentration assumption, k~ 2 0(k p )~ 1 (mf) is typically a 
factor fc 1 / 2 smaller than the leading terms. The normalization constant c\ is also 
irrelevant to the order of accuracy we seek (we will instead normalize numerically). 
Thus, keeping the leading terms in ([7^]). 

/ « ±(D 2 - D - (x ■ nf(A dn + kl))f . (76) 

At this order we also need to consider linear variation in /, so we approximate 
f(k p ) by substituting (f75j) into the / formula in ([71|) . that is, 



A second-order / expansion about k p then gives f p = f + £ P f P — (ep/2)/, which we 
can simplify]! noting the sign change in the D 2 term, to the improved 

~2 

"quadratic / estimator": /„ := f+^Df+-^(D 2 +D+(x-n) 2 (A dn +kl))f . 

(77) 

As before, we use the best available e p estimate. In d = 2 we approximate A^n/ = 
duf via a spectral differentiation matrix. Finally we normalize f p numerically. 

Figure [6] (data shown by crosses) shows the improved accuracy of this estimator: 
it gives typically one extra digit over (|75|) . and up to four extra digits over the 
trivial estimator. This error data is also summarized in the last two columns of 
Table[TJ The figures strongly suggest an empirical error of 0(e 2 fc -1 / 2 + c'k 1 / 2 ) for 
this estimator. As expected from the above discussion, the first term is a factor 
fc 1 / 2 smaller than the error of (|T5[) . As with the linear eigenfrequency estimator, 
the cubic term dominates for larger frequency distances e 3> fc _1 , which are needed 
anyway in d = 2 to capture more than O(l) mode per frequency window. Thus, 
in the fast regime, this method has asymptotic eigenfunction error 0(e 3 fc 1 / 2 ). If 
it is desired to keep this error bounded as k — > oo, one must choose e < fc -1 / 6 
rather than the e = 0(1) allowed for bounded eigenfrequency error. This reduces 
the speed-up factor of the NtD method slightly from 0(N) to 0(N 5 ^ 6 ) in d = 2. 

In section [B] we will give rigorous estimates on the Riccati estimator (|65[) and 
linear / estimator (|75)l . assuming that the spectrum of 0(fc») does not concentrate 
near j3 p . 



6 Note that one may view our procedure as inverting a Taylor series to second order. 
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Figure 7. Eigenfunctions of the pentafoil domain, (a) mode 
from the two-dimensional eigenspace at kj = 300.005956478458 • • ■ 
(function is not Z? 5 -symmetric) ; note scarring on a periodic orbit, 
(b) simple (hence Z^-symmetric) mode at kj = 300.03832269 
(All digits believed correct.) Density shows \4>j\ 2 , with white being 
zero. Parameters are as in first row of Tabled 





time / mode (sec) 
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9000 


51 
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Table 2. Runtimes and errors for the pentafoil domain of Figs. [7] 
and [8] Details are as in Table HJ except that in the reference 
method the absolute kj errors were only around 10 -7 , so errors 
below this are not discernable. For degenerate pairs, subspace 
angle replaces L 2 -error (see Remark [5.3p . Dashes show experiments 
not performed, and asterisk (*) indicates estimated ratio. 



5.3. Performance in a domain with abundant degeneracies. Having con- 
structed higher-order estimators and tested them in a nonsymmetric domain, we 
now apply them to a "pentafoil" domain parametrized by r{9) = 1 + O.3cos(50). 
For group-theoretic reasons (it has the dihedral symmetry group -D5), its Dirich- 
let eigenfrequencies are generically either simple or a doubly-degenerate pair. As 
before, we found that an N of around 6.3 points per wavelength on dfl gave full 
double-precision accuracy in computing the spectrum of 0(/c*). Table [2] summa- 
rizes our experiments comparing the proposed NtD method (using the Riccati (|65|) 
and quadratic (|77|) estimators) against the reference solver of App. |B] Observe that 
error levels of the NtD method are similar to those for the previous domain, thus 
degeneracies seem to have no deleterious effect on error. 
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Figure 8. An eigenfunction of the pentafoil domain in the two- 
dimensional eigenspace with kj = 1000.00302930323 • • ■ (all digits 
believed correct). There are around 400 wavelengths across the 
domain. Parameters are as in second row of Table [2] 



Remark 5.3. For simple eigenvalues, as before, the L 2 error \\fj — fj\\ was measured, 
with \\fj\\ = \\fj\\ = 1. For p-fold degenerate eigenvalues we used its generaliza- 
tion, the principal angle between subspaces. Here, one subspace is the eigenspace 
computed by the NtD method, while the other is that computed by the reference 
method. In the small angle limit and p = 1 this is equivalent to the L 2 error. 

Note that the reference method was slower, by roughly a factor of three com- 
pared to the nonsymmetric domain of Fig. [TJb), due to the difficulty of resolv- 
ing eigenfrequency pairs. (Here a large tolerance tol = le-6 was chosen to limit 
this slow-down.) In contrast, the NtD method pays no penalty for close or ex- 
act degeneracies — this is one of its main advantages — thus its speed-up factors are 
around three times better than for the former domain at similar frequencies. 

In Figs. [7] and [8] we show some eigenfunctions coming from the calculations of 
the first and second rows of Table [2] respectively. In the latter case, forming and 
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diagonalizing the matrix of size iV = 9000 took 3.6 hrfQ and returned 51 modes. 
Based on the previous domain, we expect mediam error similar to those in the 
first row of the table. However, since for the mode shown, kj is so close to fc* 
that we expect k error to be limited by machine precision (10~ 13 absolute error), 
and eigenfunction L 2 error to be 10 -7 . We did not attempt to run a reference 
calculation here (it would have taken 3 weeks), but using the 0(N 3 ) scaling from 
our previous tests, we estimate that our method is faster than the reference method 
by a factor of 10 3 . 

To create Fig. [SJ evaluation of the representation (|27p on a grid of 8.2 x 10 5 
points took only 27 sec per eigenfunction using the Helmholtz fast multipole method 
(FMM) implementation of Gimbutas-Greengard [37] • The entire eigenmode calcu- 
lation and plot is done by the following MPSpack code: 

s = segment . smoothstar (9000, 0.3, 5); 
d = domain(s, 1); s.setbc(-l, 'D'); p = evp(d); 
o.eps = 0.1; o. modes = 1; o.khat = 'r'; o.fhat = 's'; 
p. solvespectrum( [1000 1000. 1] , 'ntd' , o) ; 

o = [] ; o . inds = 1; o.dx = 0.002; o.fmm = 1; o.col = 'bw' ; showmodes(p 
The third line selects the Riccati esimator for k and the quadratic estimator for /. 

Remark 5.4. If a domain has a known symmetry (such as the D§ symmetry of 
this pentafoil example), it is possible to reduce N by desymmetrizing and finding 
eigenfunctions in each symmetry class separately [6, 8 . This is often done in high- 
frequency studies [9] because it increases efficiency by a significant factor. For 
simplicity, we did not implement that here. 

6. Error analysis of higher-order methods 

In this section we specialize to the case of two dimensions, d = 2. This allows us 
to use the exploit the relatively large mean spacing of Dirichlet eigenfrequencies in 
two dimensions (relative to higher dimensions). 

6.1. A spectral nonconcentration assumption. All error estimates in this sec- 
tion will be conditional on the following assumption: 

Assumption 6.1 (Absence of Spectral Concentration at scale rf). Let 77 be a positive 
real number, and let j3 be a negative eigenvalue of O(fc) satisfying 

We say that there is absence of spectral concentration at /? at the scale t) if (3 is the 
only eigenvalue of 8(fc») (counted with multiplicity) in the interval 

[£- Jp rmn(/? + |r,0)]. 

This implies, in particular, that f3 is a simple eigenvalue. 

Notice that the eigenfrequencies of A are spaced ~ 1/fc apart on average when 
d = 2, so in view of (|45l) . the eigenvalues of Q(k*) are spaced ~ 1/k 2 apart on 



This resulted in some swapping of RAM to hard drive, indicating that this about the largest 
N that can be handled on this 8 GB machine. 
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average. Therefore, for sufficiently small 77, we can expect that typically Assump- 
tion l6.1l is satisfied for most eigenvalues j3 of 0(/c*) in the range [—e/k, 0], uniformly 
in k. We will always assume that rj < 1 in our estimates below. 

One simple consequence of (|52[) is that, if e is not too large relative to 77, As- 
sumption [6TT| implies that the eigenvalue branch /3 p (k) is well-separated from neigh- 
bouring branches on the whole interval [fc*, k p ], where f3(k p ) = 0: 

Lemma 6.2. Assume that d = 2 and that the eigenvalue f3 p of 0(fc*) lies in the 
interval [— e/k, 0] and satisfies Assumption \6. 1\ at scale 77. If rj satisfies 

V > 8Ce 2 , 77 > 8Ce 3 fc, (78) 

then /3 p (fc) satisfies Assumvtion \6.l\ at the scale rj/2 for all k € [fc*,fc p ]. Here the 
constant C is the implied constant in (j52|) . 



Proof. Let f3 p (k) and f3 q (k) be two negative eigenvalue branches of O(fc), with 
f3 q (k) > ft p (k). Then (disregarding the trivial case in which /3 q (k») = 0) we have 

|&(fc*)-Ar(fc*)l > 

We will show that 

\/3 P (k)~(3 q (k)\>^L 

for all k > k* for which both eigenbranches are defined (i.e. such that both f3 p (k) < 
and f3 q (k) < 0). Using (52) we have 

d ;08p(fc)-&(fc)) 



dfc' 

since |/3 p (fc)|, |/3g(fc)| < e/k for all k > k*. Integrating over the interval [fc*, k] which 
is no bigger than e for k < k p , we find that 



e 2 e 3 , 



\[3 p (k) - [3 q (k)\ > \[3 P (K) - f3 q (K)\ - 2C VA , , , 
~ 2fc 2 ' 

using the conditions (|78|) in the last step. □ 

6.2. Error estimate for second-order eigenfunction reconstruction. We 

now give an error estimate for the estimator (1751) . given Assumption 16. II at scale 77. 
We begin with 

Lemma 6.3. Assume d = 2. For ko £ [fc^.fc* + e] and —e/k < /?(&*) < 0, the 
k-derivative of f satisfies 

f = \(W + m)f + oC-±^) (79) 
k 77 

where W is as in (I47p and to as in (|61[) . 

Proof. Consider the terms on the right hand side of (|73|) . Indeed, using Lemma l4~4l 
(and Remark 14.51 for the higher order derivatives), and since \/3\ < e/k, we see that 

< J, ^\\{x-nf{^ dn + k 2 )f\\ <e. 
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Next, using Proposition ID.2I and Assumption 16.11 at scale rj, we can estimate the 
remaining terms on the right hand side of (|73|) as follows: 

!< 



|(0 - P)- 1 (p(x ■ n) 2 (A 9n + k 2 ).f - PW'S - mf) 



< C^—p(x ■ nf(A dn + k 2 )f - pw'f ~ m f\\ LHa n) 
+ C^\\(3(x ■ n) 2 (A an + k 2 )f - pW'f - mf\\ Hl{9Q) 

= ^(\(3\k 2 + m + 1) + ^(i/3|fc 3 + \m 2 + k) 

< - e 
rj 



e 2 k) 



Here we used Lemma WM and Remark 14.51 to estimate the L 2 and H 1 norms in the 
second and third lines. Finally, the term cf is the result of projecting orthogonally 
onto the subspace orthogonal to /, so this term does not increase the norm. We 
conclude ([75]). □ 

This leads to 

Proposition 6.4. Suppose that d — 2 and that /3(fc«) = /3 p (fc*) satisfies Assump- 
tion \6.1\ at scale rj. Let f p = (x ■ n)d n (f> p /\\(x ■ n)d n <t> p \\. Then the estimator (|75|) 
for f{k p ) satisfies 



\\fp - IpWh^idn) — o(- 



(80) 



Proof. To do this, we consider two flows. One is the eigenfunction flow (|73|) above. 
The second, for the function g — g(k), is the linear flow starting at g(k t ) — /(&*), 
and flowing according to 

g = \{Wf{K)+mf{K)). 
k 

Note that the RHS here is independent of k (apart from the 1/fc prefactor). Now 
we estimate the difference between g(k) — g(fc*/(l + /?)) and the weighted normal 
derivative of the corresponding Dirichlet eigenfunction. This is a sum of two terms: 
one arising from the difference between g(k) and f(k), and one arising from the 
difference between f(k) and f(k p ), where k p is the true eigenvalue. Using (|46"j) and 
Lemma HOI together with Lemma l4~4l to see that ||W/||/fe = 0(1), the second error 
term is 

0(- + e 3 x 01 + = O - + e 3 + — + - + , 

k rj k krj rj rj 

which is certainly bounded by (|80|) for large k and r] < 1. 

The first difference can be estimated as follows. We have, with /* = /(/c*), 

^(/- 9 )4((v>-- + .»)(/-/.)) + o(i^) 

= I((W + m)(/ - 9)) + ^((W + m)( 9 - /,)) + 0(i±^). 
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Since g — /* = 0(e) (and using Lemma l4~4l again) the second term can be absorbed 
in the error term, and we get 

Again applying Lemma |4~4I we have (W + m)(f — g)\\ < C\\f — g\\. Therefore, 

K e " C(fc " fep) ll/-fflU 2 ) = 0( € -±^)e- c ^l 
Since (/ — g){k*) = 0, this inequality integrates to 

\\{f - 9){k)\\ L * =0( e —^) for \k-K\ < e. (81) 

In particular, 

||(/- 5 )(fc p )|| i2 =0(f!±i^). (82) 

V 

a 

6.3. Error estimate for the Riccati eigenfrequency estimator. Here we de- 
rive an error estimate for the higher-order eigenfrequency estimator of section 15.11 
given Assumption 16.11 at scale rj. 

Proposition 6.5. Let the frozen frequency be k z — h*. Then the estimator (|65p 
for the eigenfrequency k p satisfies 




Note that if we work in a regime with e = 0(1), then in the high frequency limit 
the term 0(ke 6 /rj) dominates in this estimate. However, as we showed in section lSTTl 
empirically the dominant error is only 0(e 5 ). When k z = |(1 + (1 + /3*) _1 )fc* is 
used instead of k z — empirically the 0(e 4 /fc) term is also absent, reducing errors 
slightly at intermediate e values. 

Proof. Consider the error in estimating the right hand side of (f62|) by (l64l) . We 
compute 

±\\(x ■ n)f\\ 2 = 2{f, (x ■ nff) = \{{W + m)f, (x ■ nf f) + o(-±^). (84) 

By integrating by parts, we see that 

1 
k 

(uniformly in e and k). Therefore, for any ko £ fc* + e], the difference between 
the value of k/3 2 \\(x ■ n)f\\ 2 (i.e. the second term of (|6"2"1) ) at ko compared to the 
value at is bounded by 

A similar calculation shows that the difference between the third term of (|62|) at ko 
compared to the value at fc* is again 0{e 3 /k + € 4 /{kr)) + e 5 /i]). Treating the fourth 
and last term of (|6"2"j) similarly, we obtain an error estimate of 0(e 2 /k 2 + e 3 / (k 2 rj) + 
£ 4 /(kn)) between the value of this term at ko compared to the value at fc*. 



((W + m)f, (x-n) 2 f) 



C 
< — 
~ k 
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Therefore, the error term in j3 for any ko £ [A:*, fc* + e] is bounded by integrating 
this error on the interval [fc*, £;*+e] and is therefore bounded by 0(e 3 /k 3 +e 4 / (k 2 rj) + 
e 5 /(kr/) + e 6 /r]). Finally, since the derivative d/3/dk is comparable to 1/k (say, 
between 1/2/c and 2/fc) by (|52p . we see that the error in the estimate for k p is 
bounded by (83]). □ 

7. Connection to the scaling method of Vergini-Saraceno 

Our above proposed NtD method is closely related to, indeed inspired by, the 
scaling method of Vergini-Saraceno [61] . Here we explain briefly the latter, using 
the language of numerical mathematics (the original paper is very short and written 
in a physics style), thus improving upon previous understandings jSJ|5]. We at least 
hcuristically explain its observed accuracy, and highlight the many differences with 
the present NtD method. 

7.1. Sketch of the scaling method. The method exploits the fact that a scaled, 
or dilated, Hclmholtz solution is still Helmholtz. Let §(k), k ^ kj, be the oper- 
ator mapping Dirichlct data to the dilational derivative of its interior Hclmholtz 
extension, that is, given g £ H 1 (dCl), and u satisfying its Dirichlet problem ([3]), its 
action is 

$(% = x ■ Vu\ dQ . 

Now consider a Dirichlct eigenfunction cf>j, and let f — (x ■ n)d n (f>j\gn- Take a 
frequency fc* — kj — e where e > is small, and define cf>j(k^) to be the dilation of 
the function <f>j to this new frequency that is 

<pj(k*){x) :=<pjfex) , x£R d (85) 
kj 

Then to first order in e, we have 

4>j(M\aa = ~0- + O(e))f (86) 

and 

x- Vc(>j(K)\ aQ = (l + 0(e))f . (87) 

The last two equations tell us that the dilated eigenmode (f>j{k*)\dn, is an approxi- 
mate eigenfunction of $(&*) with approximate eigenvalue —k/e. 

The scaling method uses a linearized self-adjoint version of the above. Let 
be the adjoint of $ with respect to (H]). The eigenvalue problem used is (analogous 
to ©), 

— ($(fc*)+$*(fc*)U = /i/i . (88) 

Although not stated as such, this is solved with the Galerkin method [37l Sec. 13.5] 
using a set of (MPS) global basis functions : CI — > C, i = 1,...,N, each 

satisfying (A + fc 2 )^(fc) = in CI. The original basis choice was plane waves 
(which seem to require CI to be convex |S]); since then, fundamental solutions [9 
and Fourier-Bessel wedge solutions [11] have also been used to handle nonconvex 
domains with one singular corner. The action of $ on ^|an is known analytically 
because each ^ is an interior Hclmholtz solution. Then the Galerkin approximation 
to ([88]) is the generalized eigenvalue problem 

GhM = ^ N) Fh^ , (89) 
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where the 'mass' matrix F has elements F,j = (£<,£j), and G has elements Gy = 
({£i,x ■ V^j) + (x ■ V£i,£j))/fc*. Further assuming that £i(k)(x) = ^i(kx), x e ft 
i.e. the basis /c-dependence is dilational, one may then check that G = dX? / dk\k=k, j 
explaining Eq. (2) of [61]. In practice, it is well known that good global bases are 
highly ill-conditioned PS] [42], thus F and G share a numerical nullspace. Then 
(|89|) must be regularized, e.g. by projection onto the numerical range of one of the 
matrices, in a similar fashion to [15 ] HO j . 

Reconstruction of eigenfrequencies is via k = — 2//i, and empirically has 
accuracy 0(e 3 ) [5], not the 0(e 4 ) claimed in [ST]. Eigenmodes are reconstructed 
from the corresponding eigenvector components fti by "undoing" the dilation via 



X hi boundary error 



U 2 (9fi) is then dominated empirically by 



0(e 3 ) with unknown fc-dependence J5[ sec. 6.3]. 

7.2. Connecting scaling and NtD methods via dilation. In place of ([87]) one 

could instead write 

(x-n)d n ^(K) = (l + 0(e))f , 

which, with ([So]) , tells us that / is an approximate eigenfunction of 0(fc») with 
eigenvalue — e/fc. It is this that motivated the authors to consider the weighted 
NtD flow — arguably more closely related to spectral theory of the Laplacian on 
n — as an alternative to dilation. 

To connect the eigenfrequency estimators of the methods, we note that $(fc) = 
O(fc) -1 + W, where W is the tangential vector field in (14"7|) . and hence that $(fc) + 
$(&)* = 26(fc)~ 1 - to, where to is defined by (1611) . Thus the operator appearing in 
([88]) can be written as (2/k*)Q(k t ,)~ 1 plus k~ x times a multiplication operator; this 
shows that the eigenvalues of ([88]) and 0(fc*) are related by /i = (2/fc*)(/3 -1 +0(1)) 



as fc* — > kj . Thus one predicts that the scaling method has eigenfrequency accuracy 
no better than that of ([2D]) ; this is observed numerically. 

For eigenfunction error, the authors are not aware of an explanation of why in 
the scaling method the combination of (|88[) and reconstruction by dilation has error 
as high-order as 0(e 3 ), as opposed to the naive O(e). Presumably the spectral flow 
of ([88]) is very close to the flow with k of (f>j(k)\gn under exact dilation. However, 
we may also connect our quadratic NtD estimator (177]) to this exact dilational flow. 
Let u be a Helmholtz solution, and let / := (x • n)u n and g :— it Ion be Cauchy 
data for its dilation u(k) to frequency k. Then one can check that / and g satisfy 
a second-order evolution equation on dfl of the form 

d 



dk 



= L 



f 



where 

r= \Ln L 12 ] = }_\ W 1 

[L21 L 22 \ k[-(x-n) 2 {A dn + k 2 ) + W W + i 

From this we can derive the first and second derivatives of / when g = 0: 



m)f; 



f=-j-2(( W + m )f ~(W + m)/ - (x ■ n y(A dn + k*)f 



(90) 



Comparing to (|74[) , we can see that the first derivative for this dilation flow at (3 = 
agrees with the first derivative for the 0(fc) flow, up to an irrelevant normalization 
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term. Moreover, the second derivative terms agree to highest order (if we agree 
that the _1 (m/) term is lower order as per Remark ID.3|) . Consequently ([77]) 
corresponds to the dilation flow just as well as it does for the O(fc) flow. 

7.3. Advantages of NtD method over the scaling method. Although the 
NtD and scaling methods have similar eigenfunction error, share the same O(N) 
acceleration factor and both are restricted to star-shaped domains, the NtD method 
has several advantages: 

• Higher-order accuracy in eigenfrequencies is possible (see section [5. ip . giv- 
ing 3 to 5 extra correct digits in practice (Fig. [5]). 

• Modes are reconstructed via ([77]) . without recourse to dilation (the latter 
requires continuation of basis functions to a strip lying outside of Q). 

• A formulation in terms of the NtD operator allows rigorous estimates such 
as Propositions 14.21 16.41 and 16.51 

• The NtD method, as implemented in Sec. 02 is robust at all choices of k*, 
whereas the scaling method is known to lose accuracy as £;* approaches 
each Dirichlet eigenfrequency [611 [5] . 

• Regularization of the numerically-singular pencil (|89[) requires a choice of 
small parameter that is not fully understood [61] HI El [IT] . 

• Our method leverages known spectrally-accurate discretizations of bound- 
ary integral operators, whereas the Galerkin method ()89[) implicit in the 
scaling method is limited by the accuracy of an available global MPS 
Hclmholtz basis. Success of the latter basis is ad hoc and quite particu- 
lar to the shape of fl. 

However, on the last point, we note that some global bases are much more efficient 
than BIE because they need only 2-3 degrees of freedom per wavelength on the 
boundary [6~Tl 19], and can be faster to evaluate than Hankel kernels. 

8. Conclusions 

We have presented, analyzed, and tested a fast algorithm for computing high- 
frequency Dirichlet eigenvalues and eigenmodes of smooth star-shaped domains in 
Mr. The acceleration is achieved by linearizing, over a frequency distance e, the flow 
of the spectrum of the weighted NtD map. The choice of weight function (x-n) -1 is 
crucial since it equalizes the gradients in this flow and prevents "avoided crossings" . 
e controls both the total time to compute all modes lying in a given frequency in- 
terval, and their resulting errors. Windows of size e are handled independently; the 
scheme is "embarrassingly parallel" . Maintaining bounded absolute eigenfrequency 
errors, one may choose e = O(l), giving a speed-up of 0(k d ^ 1 ) = O(N) over stan- 
dard methods, and more robustness since no root-search is needed. This factor is 
in practice in d — 2 roughly the number of wavelengths across the domain; we show 
an example where it is 10 3 . 

We proved robustness (neither spurious nor missing modes, see Remark 14. 7p . 
a leading third-order absolute accuracy in eigenfrequencies, and, given a spectral 
nonconcentration assumption, third-order L 2 -errors of mode boundary functions. 
This required developing some new results in the analysis of elliptic PDE of interest 
in their own right. Understanding the NtD spectral flow led to improved estima- 
tors that are empirically fifth-order for eigenfrequencies, and third-order for modes 
(with constant improved by factor k 1 / 2 ). Our scheme has many advantages over 
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the scaling method (see section f7.3j) . including an integral operator formulation, 
rigorous error analysis, and much smaller eigenfrequency errors. 

It is important to realize that the acceleration mechanism works at the operator 
level, and is therefore independent of any further acceleration that could be applied, 
such as: block iterative solvers to extract the small negative matrix eigenvalues (we 
used exclusively dense direct solvers in this work), and fast multipole or fast direct 
solvers to apply or compress the discretized operators. However, since we are in a 
high-frequency regime (oscillatory kernel), it is not at all obvious that fast solvers 
will make much difference; testing this is an obvious next step. 

Other natural questions for future work include the following: 

• Can the method be modified to remove the star-shaped restriction? 

• Can a modified method (possibly using ideas from [13]) handle other ho- 
mogeneous boundary conditions such as Neumann and Robin? 

• What accuracy can be reached for domains with corners in d = 2 or d = 3 
using appropriate BIE discretizations? (Note that the scaling method has 
been used with nonsmooth boundaries [5T1 M ITT].) 

• Can boundary error bounds on / be extended to <jp. (see Remark [32]). 

• Can ([77]) be analyzed, or improved upon in practice, while preserving the 
0{N) speed-up? One idea along these lines is high-order extrapolation from 
a e-grid of k* values; analysis would need the spectral flow for complex k. 

The reader is encouraged to try out the algorithms presented here by download- 



ing MPSpack from http : //code . google . com/p/mpspack 
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Appendix A. Smoothness of eigenvalues and eigenprojections in k 

We are interested in the flow of eigenvalues and eigenprojections of the opera- 
tor O in the parameter k. The operator O has a pole whenever k 2 is a Neumann 
eigenvalue of ft, and we wish to show that small negative eigenvalues and eigenpro- 
jections flow smoothly across such values of k. To do this we consider the Cayley 
transform of O, as in |7J). Recalling (|35l) in the case rj = 1, and solving for / and g 
in terms of u n and u\qq, we see that R(k)f = g is equivalent to the existence of u 
such that 

(A + fc 2 )u = 0, f = iu\ a n - (x ■ n)d n u\ea, g = iu\ dn + (x ■ n)d n u\ d n- (91) 

Proposition A.l. There is a neighbourhood U C C of the positive real axis such 
that there is a unique solution to the problem 

(A + k 2 )u = in Q, iu\on - (x ■ n)d n u\on = f (92) 

for every k £ U and every f g L 2 (dfl). Moreover, the solution u — u{k) depends 
holomorphically on k for k € U. 
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Corollary A. 2. The Cayley transform R(k) o/0(fc) is analytic in a neighbourhood 
U C C of the positive real axis. 

Before we give the proof of this proposition we need a couple of preparatory 
lemmas. 

Lemma A. 3. There is a neighbourhood U C C of the positive real axis such that 
for k € U , the equation 

(A + k 2 )u = in ft, iu\ an - (x ■ n)d n u\ a n = (93) 

has only the trivial solution. 

Proof. Write k — a + ib with a, b real. If u satisfies (f9"3"| . then we have 



-k 2 [ \u\ 2 + [ |Vu| 2 = [ {Au)u + [ Vu-Vi 



(94) 



1 1 1 2 



d n u u — / i(x ■ n) \u 
dfi Jen 



Taking the imaginary part we find that 



2ab I \u\ z = / (x-n)- L \u\ z . (95) 



On the other hand, we can express u in f2 via Green's representation formula (|25|) . 
It is standard that S{k) and V{k) are bounded operators from L 2 {dVL) to L 2 (il), and 
it is straightforward to check that their norms are uniform in k on compact subsets 
of the fc-axis. Using the boundary condition for u to replace d n u by i(x ■ n)~~ 1 u in 
([25]) . we see that this gives 

INfc)IU 2 (fi) < C{k)\\u\\ L 2 {dn) 

where C(k) is uniform on compact subsets. But if we combine this with (|95|) . 
then we see that for \b\ small enough compared to a, then (|9"5")) has only the trivial 
solution u = 0. □ 

Lemma A. 4. There is a compact operator L : L 2 (il) — > L 2 (il) such that Lz is the 
unique solution u to the equation 

Au = z in £1, iu\gn — (x ■ n)d n u\dn = 0. 

Proof. We define operator L\ to be inverse operator to the Dirichlet Laplacian on 
f2, i.e. L\z is the function u\ such that Au\ = z in Q with Ui|gfj = 0. It is standard 
that L\ is well-defined and compact. We then try to solve 

Au 2 = in Q, iu 2 \dn - (x ■ n)d n u 2 \an = (x ■ n)d n ui\gn; (96) 

then u\ + u 2 is the solution u that we seek. Notice that (EMI) implies that 

(x ■ n)<9„ui|ao = (i- B(0))u 2 \dn, 

where B(0) = (x ■ n)A(O) is the weighted Dirichlet to Neumann operator at zero en- 
ergy. The operator B(0) is self-adjoint on L 2 (dfl) with our weighted inner product, 
so we can invert i — B(0) and find that 

u 2 \on = (i~ •S(0)) _1 ((a; • ri)<9 n wi| a n). 

Finally, if P is the classical Poisson operator taking functions on dft to the harmonic 
function in f2 with the given boundary value, then we have 

u 2 = P o{i- B^S)Y Y {x-n)d n L 1 z. 
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We recall some standard mapping properties of these operators. The operator 
L\ maps L 2 (f2) to iJ 2 (f2), then (x ■ n) times the normal derivative of this at the 
boundary maps to i/ 1 / 2 (90), then (i — _B(0)) _1 is a pseudodifferential operator 
of order -1, hence maps H 1 / 2 ^) to H 3 / 2 (n), while P maps H Z / 2 (VL) to H 2 (ft) 
[53l Ch. 5, Prop. 1.7]. Denote the composite operator L2, i.e. 112 — L2Z. Then 
we see that L2 maps L 2 (f2) continuously to H 2 (Vl) 1 and hence, using the compact 
embedding of H 2 (fl) into L 2 (f2), we see that L2 is compact on L 2 (il). Hence 
L = L\ + L2 is compact. Uniqueness of the solution follows from Lemma [A. 31 with 
k = 0. This completes the proof of Lemma [A. 41 □ 

Proof of Proposition. Let U be as in Lemma lA. 31 Then this lemma guarantees the 
uniqueness of u satisfying (|92j) . It remains to establish existence. To do this, we 
first find W\ such that 

Atui = in VI, iw\dn - (x • n)d n wi\an = f 

which is done exactly as in (f96|) . Then we look for W2 satisfying 

(A + k 2 )w 2 = -k 2 wi, iw 2 - (x ■ n)d n w 2 \dn = 0. 

If we can find such a W2, then u = Wi + W2 is our required solution of (]92[) . Using 
the operator L from Lemma lA. 41 this can be written 

w 2 = -L(k 2 w 1 + k 2 W2), 

which is equivalent to 

(ld+fc 2 L)w 2 = -k 2 Lwi. 

Thus we get a solution provided that Id+fc 2 L is invertible. Since L is compact, this 
will be the case provided that Id +k 2 L has trivial null space. But if v is in the null 
space of this operator then v satisfies (|93|). which means by Lemma [A. 3 1 that indeed 
v = 0. Therefore the null space is trivial, so Id+fc 2 L is invertible and we can find 
W2 as above. This establishes existence of u. Finally, using the compactness of L 
and analytic Frcdhom theory Thm. VI. 14], for k € U, (Id +k 2 L)- 1 (-k 2 Lw 1 ) 
is analytic in k, showing that u(k) is analytic in k. □ 

It follows from the analyticity of R(k) that in any interval I of the unit circle 
in which the spectrum of R(k) is discrete at k = k , the eigenvalues of R(k) in 
I are analytic as a function of fc, and one can choose an orthonormal basis of 
the corresponding eigenspaces that varies analytically 34, Ch. VII, sec. 3]. This 
implies that the eigenspaces of O(fc) vary analytically in any interval where the 
spectrum is discrete, with the exception of a finite number that have a pole at each 
Neumann eigenfrequency. Since Q(k) is a pseudodifferential operator of order —1 
and therefore compact, this means that the eigenspaces vary analytically except 
when the eigenvalue hits zero. In fact, we can say more. Before we state the next 
proposition, recall that the eigenvalues of 0(fc) are monotonic increasing in k — 
see ([17]). 

Proposition A. 5. (i) Suppose that is an eigenvalue o/0(fc*). Then there is an 
analytic eigenvalue branch /3(fc) with (3 f as k f k*, and the multiplicity of the 
eigenspace is equal to the sums of the multiplicities of all such branches. 

(ii) Conversely, suppose that P{k) is an eigenvalue branch of<d(k) tending to zero 
as k t fc*. Then k* is a Dirichlet eigenfrequency, the eigenprojection has a limit 
as k t k*, and it is the projection onto a subspace of the space of weighted normal 
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derivatives of Dirichlet eigen] 'unctions with eigenfrequency fc*. The eigenvalue P(k) 
is C 1 as a function of k up to and including k — fc* , and satisfies (I19|) . Finally, 
if the eigenvalue /3(fc) is simple up to and including fc = fc*, then the eigenfunction 
/(fc) is C 2 up to an including k — fc*, and satisfies (I74|) . 

Proof, (i) Suppose that is an eigenvalue of 0(fc*), with eigenspace V. Choose 
an interval (a, b) containing 0, with neither a nor b in the spectrum of 0(fc*), and 
such that there are no eigenvalues in the interval (a,0), and let II a .fc(fc) denote 
the projection onto the eigenspaces of 0(fc) with eigenvalues in the interval (a, b). 
Define 9(fc) = II a .fc(fc)0(fc). Then for fc close to fc*, 9(fc) is an analytic family, 
again using [MJ Ch. VII, sec. 3]. By the calculation in Lemma \2.2l dQ(k)/dk is a 
positive operator. Since 

(e(fc*)/,/) =0foraU/e V, 

we have 

(&(k)f, /) < for all / G V, fc < fc*. 

So there are at least (dim]/) negative eigenvalues of O(fc), which tend to as 
fc t k*- These branches are analytic for fc < fc* since the negative spectrum of 9(fc) 
is discrete. The statement that there are exactly (dimV 7 ) negative eigenvalues of 
0(fc) which tend to as fc f fc* follows from the proof of (ii) below. 

(ii) For simplicity, we first prove (ii) assuming that /3(fc) is simple. In that case, 
taking the eigenfunction /(fc) to be normalized in L 2 (dQ), we see from the identity 
(|110p that the extended eigenfunction u(k) is uniformly bounded in i/ 1 (il) as fc t 
Therefore, there is a sequence ki tending upward to fc* such that u(ki) has a weak 
limit v in H X (Q), and therefore, a strong limit in L 2 , along this sequence. It also 
follows from (fT"T)) and ((5^)) that the L 2 norm of u(ki) does not tend to zero along 
this sequence, so v is nonzero. From the fact that the u{k{) are Helmholtz, we find 
that along this sequence, we have 

lim f u(ki)(A + fc 4 2 )V = for all V € C£°(fi), 
Jn 

implying that v is a weak solution of the equation (A + k 2 )v — 0. By elliptic 
regularity, this means that v is smooth in the interior of fi and satisfies the equation 
in the strong sense. Also, using the continuous map from _ff 1 (J7) — > L 2 (dft) given 
by restriction to the boundary, we see that v\ga is the weak limit (in L 2 {d£l)) of 
u(ki)\on- But u(ki)\dn tends strongly to zero (since u(k l )\g^ = j3(ki)f{ki) and 
tends to zero) and a fortiori weakly, so v is zero at the boundary. It follows that v 
is a Dirichlet eigenfunction. We see that u(k) has a continuous extension to fc — fc*, 
such that it is a Dirichlet eigenfunction at fc = fc*. That is, is an eigenfunction 
of 0(fc*), so the eigenvalue branch /3(fc) extends continuously to fc = fc*. Given 
([51)l . we see that $ has a limit 1/fc* as fc f and therefore, /3(fc) is C 1 up to an 
including fc = fc*, and (|19p holds. 

If (3 is a multiple eigenvalue, we proceed similarly. We take a sequence of extended 
eigenfunctions as before, and produce a Dirichlet eigenfunction v\. Next we take 
another sequence of extended eigenfunctions orthogonal (at the same value of fc) to 
the first sequence, and produce another Dirichlet eigenfunction vi, and so on. We 
find a subspace of Dirichlet eigenfunctions at frequency fc* of dimension equal to 
that of the multiplicity of (3(k). 
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Again assuming that the eigenvalue f3(k) is simple up to and including fc — fc*, let 
77 be a positive number such that Assumption 16.11 holds in some interval [fc* — S, fc*] 
for some S > 0. Then Lemma 16.31 and Lemma 14. 4[ show that f(k) is uniformly 
bounded as fc f fc*, and hence /(fc) has a limit as k — > fc*. Now referring to (f73|) . 
using the continuity of f(k) just shown, Lemma T4. 41 and (|5Tj) to bound derivatives 
of /(fc), and Proposition ID. 21 to control the norm of (O(fc) — /3) _1 uniformly as 
k t we see that / itself is continuous up to k = fc* . Iterating once more using 
(1731 . we see that / is continuous up to k = fc*. Hence /(fc) is C 2 up to fc = fc* and 
formula ([73)) extends by continuity to k = fc* to yield ([74]) when /3 = 0. □ 

Appendix B. Computation of reference eigenfrequencies and 

EIGENMODES 

Here we describe our implementation of a standard published method for compu- 
tation of eigenpairs, which we use as a reference to assess both accuracy and speed 
of the NtD method. Recalling the definition (|30l) . we have the following standard 
result (e.g. see }411 Lemma 8.4] which applies for domains with Lipschitz boundary; 
note the opposite sign convention). 

Lemma B.l. A positive frequency k is a Dirichlet eigenfrequency if and only if 
the operator (| — D*(k)) has a non-trivial nullspace. Furthermore, its nullspace is 
precisely the space of boundary normal derivatives of solutions of (A + k 2 )u = in 
fi with homogeneous Dirichlet data on the boundary. 

Its proof uses the jump relations and the uniqueness of the exterior Helmholtz 
Neumann boundary value problem p~8j [41] . The numerical method is then, follow- 
ing Backer [6j sec. 3.3], to search along the k axis for (near) zeros of the lowest 
singular value of a matrix discretization of the operator (| — D*(k)). We use the 
Nystrom quadrature as in (|39)) -([40l); the same N as before may be used to achieve 
quadrature errors around machine precision. The cost of each minimum singular 
value evaluation is then 0(N 3 ). (We note that finding roots of the determinant is 
faster but is not able to distinguish close eigenfrequencies or handle degeneracies 
reliably [5]). 

The minimum singular value, which we call t, as a function of fc, has the form 
of a series of V-shapes with the bottom of each 'V approaching zero (e.g. see 
Fig. 8 of [IH] or Fig. 5.1 of [TD] ) . Reliably locating all such minima in a range 
of fc is not trivial, crudely speaking because close eigenfrequencies lead to small- 
scale W-shapes that are difficult to distinguish from a 'V. We make use of the 
empirical observation that the slope of t(k) appears to have an upper bound Ct of 
size 0(1) which depends only on fi, and that most of the V-shapes have this slope. 
We initially evaluate t(k) on a regular grid of spacing about 0.2 times the mean 
eigenfrequency spacing. At each local minimum on this grid we use the information 
about higher singular values to decide whether to i) perform fitting of a parabola 
to the three neighbouring samples of t 2 (fc), and iterate this fit procedure until 
convergence, or ii) recursively call the same routine on an (about 3 times) finer grid 
covering three (or more, if there are nearby small values of t) neighbouring grid 
points. We omit several details of the algorithm required for robustness. 

This has been coded into MPSpack and may be run (for instance for the example 
of section 13. ip via 

o.maxslope = 1.5; o.tol = le-12; p . solvespectrum( [90 100], 'ms', o) ; 
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where maxslope defines the value Ct, and tol the requested absolute tolerance on 
k. When Ct is chosen correctly, the algorithm finds all kj in a given k interval, 
needing around 15 evaluations per simple eigenfrequency found, and typical errors 
are 10~ 13 or less. When eigenfrequencies are degenerate, many more recursions 
are needed to establish reliably that they are not distinct; for instance at o.tol = 
le-6 it still requires around 50 evaluations per multiple eigenfrequency found (this 
scales like log tol), and typical errors are 10~ 7 . 

Once accurate eigenfrequencies have been found, modes are found as follows. For 
each kj , the last right singular vector of the above matrix is computed at a cost of 
0(N 3 ); according to Lemma \B. II this approximates d n <f>j at the quadrature nodes. 
Normalization is done via (p~8]> . Eigenfunctions <fi may then be reconstructed via 
([25)) . In the case of an p-fold degeneracy, the last p right singular vectors are used. 
The whole method thus scales as 0(N 3 ) per mode with a rather large constant. 

Appendix C. Proof of Lemma l4~4l 

Proof of Lemma \4-4\ We prove the theorem under very slightly more general con- 
ditions. That is, we replace (x ■ n) — both in the boundary condition (1491) and 
in the weight factor in the inner product on L 2 (dil) — by an arbitrary smooth 
positive weight, which we denote m. First, we introduce a spectral cutoff. Since we 
are using a weighted inner product we define the operator 

AaO jUJ = V^n Vtan 

on L 2 (dfl), where *' w denotes the adjoint with respect to the weighted inner prod- 
uct. (Below, we write * instead of *'"' but all adjoints in this appendix should be 
understood to be with respect to the weighted inner product.) We write 

Id = V{A dn , w /k 2 ) + (1 - *)(A ao ,™/fc 2 ), on L 2 {dn), 

where is 1 for t > 3/2 and for t < 5/4. Let us write * for ^(A dn , w /k 2 ) 
below; note that ^ is a semiclassical pseudodiffcrential operator of order (0,0) |3 
supported where \r}\ > 3/2. We can expect that the L 2 norm of (Id — ^)V t anU is 
bounded by 2k times that of u, since applying (Id — <]/) removes frequencies of order 
> 2k. To verify this, given a vector field W of unit length and tangential to the 
boundary, we compute 

\\(ld-^)Wu\\ 2 L2{9n) = \\W(ld-*)u + [1 - *,W]u||i a(8n) 

< ^\\W(ld-y)u\\ 2 LHdn) + 4|| [1 - W]u\\ 2 LHdSl) (97) 

= ~{{ld-y)W*W{l& -¥)u, u) + 4|| [1 - % W}u\\ 2 L2{my 

®A operator with parameter h is a semiclassical pseudodiffercntial operator of order (I, m) on 
Of! if its Schwartz kernel can be written locally (that is, with respect to some local coordinate 
patch y = (yi , . . . , yd—i)) in the form 

h . (d . 1) . l r e^-^^aiy^Mdr), 
where r) S R d_1 and the symbol a is smooth in h and satisfies symbol estimates 

\8^a(y, V ,h)\ < C a , 7 (0TtoP) m - 171 . 
Here the parameter h is k~ 1 . 
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Notice that [1 — 4 r , W) is a semiclassical pseudodifferential operator of order (0, — oo), 
hence with uniformly bounded (in k) L 2 (dQ,) — > L 2 (<9f2) operator norm. If we sum 
over an orthonormal basis W\, . . . , W n -i, then using J2i W*Wi = Agn,w and the 
fact that 1 — ^(t) vanishes when t > 3/2, we have 

||(M-¥)Vtanu||£ a(8 n) < 2<(M-*)A 8n ,t (Id -*)«,«> +2^ ||[1- WiHli a(en) 

i 

< (2k 2 + C)\\u\\ 2 L2{dny 

(98) 

Next we analyze the high energy part, ^Vtan - ". We use the single and double 
layer boundary integral operators S(k) defined by (|29]l. and D(k) defined by ([30)1, 
We also write Q(k) for the (hypersingular) operator d n:c d ny Go(k; x, y) restricted to 
the boundary in both variables. 

We now quote results from Section 4]. Here it is shown that S(k) and D(k) 
are pseudodifferential operators of order (—1,-1) in the 'elliptic region' {\rj\ > 1} 
(where \r]\ is the length of 77 with respect to the induced boundary metric on dSY), 
in the sense that if $ is a semiclassical pseudodifferential operator of order (l,m), 
microsupported in the elliptic region, then §S(k) and $D(fc) are semiclassical pseu- 
dodifferential operators of orders (I — l,m — 1). Moreover, the analysis from [5^1 
Section 4] applies to $Q(k) which shows that $Q(k) is a semiclassical pseudodif- 
ferential operator of order (I + 1, m + 1), with principal symbol — \ko{<&)^J\ — |??| 2 
where <t($) is the principal symbol of $. (See Remark IC. II in case this is confusing.) 

For any Hclmholtz solution u we have the Green's representation formula (|25[). 
By differentiating normally at the boundary dtt, we obtain, using (|33[) . 

d n u(x) = (D(jfc)* + \)d n u - Q(k)u. (99) 

Let us write £>{k) for the kernel m~ 1 D(k)m. We then obtain from (|99| and the 
boundary condition (|49| that 

^j3md n u = 0(D(k)) t (md n u)-l3mQ(k)u =^> u = 2(D(k)) t u-2(3mQ(k)u. (100) 

Next we differentiate tangentially, apply ^P, and take the inner product with 
^Wu, where W is a tangential vector field of unit length. We obtain 

(tfWu, ^Wu) = 2($ 2 W(D(k)) t u, Wu^ - 2/3(W*y 2 WmQ(k)u, u). (101) 

Using the results of [29 mentioned above, we see that 1 $ 2 W(D(k)) t is a semiclassical 
operator of order (0, 0), hence bounded on L 2 (dQ) uniformly in k. (Here we use the 
property of ^ that it is microsupported in the elliptic region, in fact in the region 
{\f)\ > 4/3}.) Hence the first term in (|101j) is estimated by 

C\\u\\ L 2 {dn) \\V tari u\\ L 2 {dn) . (102) 

In the second term, we have the operator W*^ WmQ(k). Since W*^ 2 Wm is 
also a pseudodifferential operator microsupported in the elliptic region, and since 
W and W* are of pseudodifferential operator (1, 1), we see that W*^ 2 WmQ(k) is 
a pseudodifferential operator of order (3,3), with principal symbol 

- fc 3 |cr(/iM/)| 2 m^ 2 (?7)Vl?7l 2 - !• (103) 
This is minus the square of a smooth symbol, namely 

k 3/2 m 1/2 ^(r 1 )a(ihW)(\f 1 \ 2 - 1) 1/4 . 
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The sign of (|103|) is crucial, as it will effectively allow us to discard this term, 
which would otherwise be too big to estimate. This works as follows: by the 
pseudodifferential calculus, we have W*^ 2 WmQ(k) = —B*B + A 2 , with A 2l B 
semiclassical pseudos, B of order (3/2,3/2) and A 2 of order (2,2). This term can 
therefore be expressed 

2/3\\Bu\\l-2p{u,A 2 u). 
Since A 2 is supported in {|?7| > 3/2} we can write 

A 2 = V* an A Vtan + Ax 

with Aq of order (0,0) and A\ of order (1,1). Since A\ can be chosen to be 
microsupported in the elliptic region, we have A\ = A' ■ V ta n + Aq, where A' Q and 
Aq are pseudos of order (0, 0). Thus we have 

-p(W*V 2 WmQ(k)u, u) = 2l3\\Bu\\ 2 2 -2(3^(V taD u, A V t a„M) + (^A^ 2 , u)M + A%u, u) 

This gives an estimate for the second term of (jlOlf) of the form 

P\\Bu\\l + C\/3\ (\\V tan u\\ 2 L2(m + \\u\\ 2 L2(m ) (104) 

(where we dropped the term C|/3||| , u||x,2(aQ)||Vtan' u lli, 2 (an) since it is controlled by 
the other terms on the right hand side). Combining (|102l) and (I104[) . we find that 

\\*Wu\\Z </3||Su||| + C|/3|(||V tan n||| 2(an) + ||u||^ (fln) ) + C\\u\\ L 2 idn) \\V ta , n u\\ LHd( 

< C|^l(ll V tanU||i2 (ao) + || U|| |,2(0n)) + C||w|U 2 (aO)l|VtanU|| L 2 (an) 

where we are able to discard the Bu term since /3 < 0. Combining this with (|98|) 
and using the inequality ||a + 6|| 2 < 3||a|| 2 /2 + 3||6|| 2 we find that 

\\Wu\\ 2 < (3k 2 + C)\\u\\ 2 L2{dn) + C\\u\\ L 2 (dn) \\\7 tan u\\ L 2 (dn) + C|^|||V tan u||| 2(an) , 

and summing over an orthonormal basis of W we find that 

HVtanWlla < ( 3fc2 + C )\\ u \\h(an) + C|l«IU 2 (9n)||V t anU||L2 (an) + C|/3|||V t anW|i 2 2 (0o) . 

Finally we write 

C\\u\\ L 2 {dn) \\\7 tan u\\ L 2 {dn) < AC 2 \\u\\ 2 L 2 (m) + ^l!V tan w|li2 (ao) , 

and observe that for |/3| < C/16 we can absorb the || V ta nw||| terms on the left hand 
side to deduce 

||VtanU||^< fa 2 + 4C 2 )\\u\\ 2 L 2 {m} . 

For k > K, we have 8(3fc 2 + 4C 2 )/7 < 4fc 2 and we arrive at (|5D|1. □ 

Remark C.l. Since D(k) involves one extra derivative than S(k), it might seem 
peculiar that S(k) and D(k) are both order (—1,-1) in the elliptic region. In 
fact, the distributional limit of d ny Go(k; x,y) to the boundary in both variables is 
a pseudodifferential operator of order (0,0) in the elliptic region, but the leading 
part of this operator is half the identity — supported at the diagonal — so it 
does not appear when the kernel function is restricted to the boundary in both 
variables; rather this part of the operator shows up as the \ in the jump formula 
(|3Tj) . This does not happen for Q(k), hence its order is two more than that of S(k), 
as expected. 
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Appendix D. Estimates involving (6-/3) 1 

To prepare for this operator norm estimate we first generalize an estimate from 
[TJ] from the Dirichlet boundary condition to the 'near Dirichlet' Robin boundary 
conditions, that is the boundary condition 

u\ 9n - P(x ■ n)d n u\en = 0, -- < p < 0, (105) 

where 6 is small and j3 < 0. This is a self-adjoint boundary condition and there is a 
corresponding orthonormal basis of eigenfunctions (f)j, with eigenvalues Ej — (kj) 2 . 
Then we have 

Proposition D.l. Let dQ be smooth. Then there exists S > and a constant Cn 
depending only on Q and 6 such that the operator norm 

!l ^ (x-n)d n <t>?({x-n)d n <f>P-) (106) 



L 2 {dQ.),{x-n)- 1 da 



is bounded by Cak 2 , uniformly for j3 in the range [— 8/k, 0]. 

Proof. We use the same method of proof as in [12] , but with additional work since 
we no longer have our functions vanishing at the boundary. 
Our starting point is the identity 

cj)[A + k 2 ,V}(t)= I {& + k 2 )<j)Vcj) 

Jn (107) 
cj)V(A + k 2 )(b+ / 4>d n V<f>-d n <i>V<l>. 



I an 

If we choose V to be a smooth vector field equal to (x ■ n)d n at the boundary, then 
the last term in (|107[) is ||(x • n)d n 4>\\ 2 . We can therefore express 

||(x-n)a„0|| 2 = / {x-n)\d n <P\ 2 
Jon 



(j)d n (x ■ n)d n <t> + / (A + fc 2 )#>- / c/)[A + k 2 ,V](f> - / <j>V {A + k 2 )(j) 
ion Jn Jn J 

= I + II - III - IV. 

(108) 

The first step in the proof of Proposition lD.ll is to estimate this squared L 2 norm 
when <j) is an approximate eigenfunction, that is, a function satisfying 

iy(n) = l, + • n)d n <p = at dfl, ||(A + k 2 )4>\\ L2(n) = O(k). (109) 

We claim that this implies that 

|| (x ■ n)d n (f>\\ L 2( dn) < Ck, 

which we prove by estimating terms / — IV in (|108[) by by a constant times k 2 . 
Before doing so, observe that using (|105p and Lemma 14.41 any boundary term of 
the form 

/ aifc 2 H 2 + a2|V tan </>| 2 + ka 3 (/)\X7 tan </>| + ka^d n <\> + a 5 a„0|V t an0| 
Jon 

where a, are bounded functions on dQ, not depending on (f> or k, can be estimated 
by C(5||9 n i^||| and therefore (for sufficiently small 8) can be absorbed in the left 
hand side; we will call them 'acceptable' boundary terms. 
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Consider the identity 

I |v^| 2 - / (-A0£+ / a„#. (no) 

Jo Jn Jan 

For /3 < 0, the last term is negative using the boundary condition (|105j) . implying 
(using ([Tog]) ) that 

HWH^n) = O(fc) =s- ||V0|U W = O(fc). (Ill) 

Using this we see that the term 77 on the right hand side of f|108[) is 0(k 2 ). 
Term IV can be expressed after integration by parts as 

f V(j}(A + k 2 )<j)+{divV(j))(A + k 2 )(t)- [ (x ■ n)(V ■ n)cj)(A + k 2 )cj). 
Jn Jan 

The first two terms are dealt with as above. In the third term, we expand A = 
d 2 + AgQ + (d — 1)77 <9„- Notice that the d 2 term cancels term 7 in (|108j) up to an 
acceptable boundary term. So we have to estimate the terms 

/ 4>{A au + (d- l)Hd n + k 2 )<i>. 
Jan 

The Hd n and k 2 terms are acceptable. The Aqq term is estimated by integrating 
by parts to convert the integrand to |Vtan0| 2 which is also acceptable. 

To estimate term 777 we use the fact that [A, V] is a second order operator and 
therefore of the form of a finite sum ViWi where Vi, Wi are smooth vector fields. 
We can integrate by parts modulo an acceptable boundary term and obtain 

f y^v^w^, 

Jan 

and the L 2 norm is bounded by C\\ Vtt||^ 2 (Q) = 0(k 2 ). This completes the proof 
that || (a; • n)d n (j>\\ 2 L2{dn) = 0{k 2 ). 

The second step of the proof is the same as in [JJ]. For the reader's convenience 
we repeat the argument here. We define an operator T from the range of the spectral 
projector hk— i,fe+i] (v~ that is, the vector space spanned by eigenfunctions of 
—A with eigenvalues in the range [(k — l) 2 , (k + l) 2 ], to L 2 (d£l), by 

T(f> = (x ■ n)d n (f>\en- 

Then, any such <p satisfies 

meaning that <f> is (after normalization) an approximate eigenfunction in the sense 
of (|109|) . So in the first step of the proof above, we showed that \\T(f>\\ < Ck\\<f>\\, 
or in other words that T has operator norm at most Ck. But the operator norm of 
7 is equal to that of 7*, and the operator that appears in (|106l) is precisely 77*. 
Since the operator norm of TT* is precisely the square of the operator norm of T, 
this completes the proof of the theorem. □ 

We now prove an estimate on the term involving the generalized inverse (6— /3) _1 
in (3). 
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Proposition D.2. Let (3, for k £ [fc*, fc* + e], be an eigenvalue of 0(fc) satisfying 
—e/k < ft < and satisfying Assumption \6.1\ at the scale r\. Assume also that 
r\ satisfies (|78[1 . Then there exists C depending only on fl such that for any z € 
H 1 (dfl), we have the estimate 

\\(Q(k) - P)- X z\\ L , [m) < c(j\\z\\ L * iaQ) + \\z\\ m{an) ). (H2) 

Proof. Let / be the eigenfunction of 6(fc) with eigenvalue f3. Let z £ L 2 (dVl) be the 
projection of z into the subspace orthogonal to /. Let v be a Helmholtz solution at 
frequency k such that (v— /3(x-n)d n v)\gn = z; then (0(fc) — (3)^ 1 z = H^(x-n)d n v\g^i, 
where Ilf is the orthogonal projection onto the subspace orthogonal to /. We need 
to estimate the norm of (x ■ n)d n v\gn relative to the norm of z; for this it suffices 
to assume z = z. 

To estimate the size of (x ■ n)d n v\gci, we expand v in eigenfunctions (p% as in 
Proposition lD.il We write — k 2 . By assumption, there is a j = j* such that 
Ej^ — E 13 ] the corresponding eigenfunction <^ satisfies (x ■ n)d n 4>^^ = f . 

Let 

We may assume that a,j w = 0, as Hjr((x ■ n)d n (f>jj = 0. Then 



Therefore, 




If we try to take the normal derivative term by term in this series and sum, unfor- 
tunately we end up with a divergent series (even when the quasi-orthogonality of 
the boundary values (x ■ n)d n 4>j is taken into account). To avoid this problem we 
write v' as the solution to 

Av' = in Q, (v - j3{x ■ n)d n v)\gn = z; 

there is a unique solution to this problem due to the negativity of (3. Then we can 
express v' = Y] a 'j<t>^ where from a similar computation to above 
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Therefore, v — v' has an expansion 



(113) 



which has improved convergence properties as the denominator is now ~ (E 



instead of {Etf) 1 , as j — > oo. From this we see that 

{x ■ n)d n {v ~v')=Yl { e P(ep_ e P) j dn zd "$Y x ■ 

Now we use Proposition ID . 1 1 proceed as in Section 4 of [12] and show that the 
operator 

has operator norm at most C+Ck 2 /d(E" , cr*), where d(E@ , a*) denotes the distance 
from E@ to the nearest point of the spectrum on A with boundary condition (|105l) . 
By Assumption 16.11 and Lemma 16.21 this distance is at least Cr], so we get an 
estimate on the operator norm of Ck 2 /rj. Therefore v — v' has norm at most 
Ck 2 /rj\\z\\ L 2 (m) . 

To treat the term v', notice that that (x-n)d n v' is (8(0) — we will estimate 
the operator norm of (6(0) — /3) _1 . The operator 0(0) is a positive operator, since 

(v, (x ■ n)d n v) = / vd n v = / vAv + \Vv\ 2 > 0. 
Jan io 

Therefore, as (3 is negative, the norm of (0(0) — (3)^ 1 z is no bigger than that of 
Q(0)~ 1 z. The operator O(0) _1 , which is nothing other than the multiplication 
opertor (a; • n) composed with the Dirichlet-to-Neumann map at energy zero, is a 
pseudodifferential operator of order 1, and therefore 

< C\\z\\ H i( dn y 

This concludes the proof of Proposition ID.2I □ 

Remark D.3. In fact, although the above analysis shows that ||(0(fe) — /3) _1 z| 
can indeed be as large as Ck 2 /d(E" , cr*) times ||z||, this only happens in a 'worst- 
case scenario' in which z is a multiple of d n (j>^ where Ej is the eigenvalue of Q(k) 
closest to (but distinct from) E@ (or, more precisely, a linear combination of an 
0(1) number of the [x ■ n)d n <p^ with closest eigenvalues). In a more 'typical- 



case scenario ', the coefficients a,j, for \Ej - E fi \ < vF, would be ~ k 1 / 2 



m 



magnitude — this can be seen from Proposition ID.ll and the arguments of [12] , 
which show that the (x ■ n)d n (j)j have norm ~ k and are approximately orthogonal 
for \E? - E p \ < yJW. On the other hand, for \E$ - E f3 \ > VE$, we gain a power 

of V E& — k in the denominator of (|113[) . This suggests that, typically, we would 
have ||(Q(fc) - P^zW no bigger than a constant times k 3 ^ 2 \\z\\/d(E^, a*). This 
would imply that in formula ([74"]) for the second derivative of / at (3 — 0, and 
given Assumption 16. ll with r\ ~ 1, the Q{k)~ l (mf) term is usually smaller by a 
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factor <~ k~ x / 2 than the principal terms, even though it is of the same order in the 
worst-case scenario. This is a heuristic justification for dropping this term in the 
quadratic estimator ([77)1 . 
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